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Synopsis 


Application of linear elastic fracture mechanics to design of machine or structural com- 
ponents is restricted by the condition that the plastic zone at the crack tip be small 
compared to the size of the K dominant field and any relevant geometric dimension. It 
is virtually impossible to satisfy this condition for high toughness, low strength materials 
which generally undergo extensive plastic deformation and crack-tip blunting prior to 
the beginning of crack growth. Crack initiation in these materials is usually followed by 
a stable crack growth before it fails in catastrophic majnner. The need to incorporate 
the influence of significant plastic deformation in crack growth initiation ( i.e. the onset 
of growth of a existing crack ) and subsequent stable growth is the main impetus for 
the development of elastic-plastic fracture mechanics. While the development in this 
field is ongoing, the theory of ductile fracture is not yet fully established. The use of 
HRR ( Hutchinson [1] , Rice and Rosengren [2] ) stress field to predict the stress and 
displacement fields ahead of a crack tip is restricted by several limitations like propor- 
tional loading, small crack-tip geometry change etc. Therefore, the use of Jjc as a ductile 
fraeture parameter becomes dependent on several conditions. 

Recent investigations on ductile fracture show that the crack growth initiation in 
a ductile material is a microscopic phenomenon occurring over a small zone around the 
craxd;-tip. This zone is named aa fracture process zone and is of the size of critical 
crack-tip opening displacement ( CTOD^ ). Within the fracture process zone, profuse 
void growth and coalescence occur which cause a local instability in the plastic flow and 
ultimately lead to ductile failure. 

Modelling of the above phenomena in the framework of continuum mechanics in- 
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voives the following aspects: 

(i) Constitutive modelling of the material behaviour to incorporate the effects of mi- 
crostructural features of void growth, 

(ii) Representation of void coalescence condition in terms of continuum paraiueters so as 
to identify a critical continuum parameter for micro crack initiation, 

(iii) Use of the above critical continuum parameter along with a length parameter ( repre- 
senting the microstructural features of the material like grain size, inter-particle spacing 
etc ) for the development of a local criterion for crack growth initiation. 


In this thesis, the problem of ductile fracture is formulated with the help of contin- 
uum damage mechanics ( CDM ). The change in constitutive relations of the material 
due to void growth has been incorporated by introducing an internal damage variable 
{D) in the plastic potential. The damage variable is physically identified as the area void 
fraction at a point. Based on Lemaitre’s CDM theory [ 3 ] and the experimental results 
of LeRoy et al [ 4 ] , a damage growth law is proposed as: 


D = + {ar + a2D){-Y)e%. 


( 0 . 1 ) 


Here, oi, 02 and c are material constants derived from LeRoy ’s experimental results. 
(— T) is the thermodynamic force corresponding to the damage variable (D). The ex- 
pression of {—Y), as derived from Lemaitre’s CDM theory, is given by 


-F = 


'eg 


2E{1 - DY 


(1 + 1-) + 3(1 - 


'eq 


( 0 . 2 ) 


The ratio (^) is called as triaxiality where cr^ is the mean stress and a^q is the von 
Mises equivalent stress. E is the Young’s modulus, u stands for poisson’s ratio and 
is the equivalent plastic strain rate. Introduction of internal damage variable in plastic 
potential modifies the plastic flow rules. The modified plastic flow rules are stated below: 


•p _ ^3 o'ij 

1-D2cr,,’ 

(0.3) 

A 

, ( 0 . 4 ) 
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Here, cr- is the deviatoric part of the stress tensor (T,j, p is the har<d<;ning variable and A 
stands for a scalar multiplier. 


The constitutive equations cited above are used in an updated Lagrangian elastic- 
plastic FE programme which includes both material and geometrif nonlinearities. The 
condition of void coalescence is derived from Thomason’s [ 5 ] limit load model with ap- 
propriate modifications to suit the FE programme. The expression of the above condition 
is given below: 


CTl 


0.1 -I- 


L2 

{1 - exp{-elqj2)}2 _ 


= 0 . 


(0.5) 


Here, ctj is the maximum principal stress and Cy stands for current yield stress of the 
intervoid matrix material. Whenever the above condition is satisfied during the defor- 
mation process, micro crack initiation is assumed to take place. 


Different case studies are performed with specimens of varying geometries. AISI- 
1090 spherodised steel is taken as the test material. It is found that, at micro crack 
initiation, the damage \-ariable (D) reaches a critical value (Dc) which is almost constant 
for all the cases. Hence, the critical value of damage variable [Dc) is taken as a material 
property independent of specimen geometry. Further, it is observed from the case studies 
that a ductile fracture process depends on both plastic strain and triaxialjty (f^)- It 
should be noted that the damage variable incorporates the role of both plastic strain and 
the triaxial ity ( equation (1) ). 

Once the critical damage (Dc) is established as an aecepted critical continuum pa- 
rameter for micro crack initiation, it is used in a local criterion, to predict the crack 
growth initiation in Mode-I ductile fracture. Average austenite ^ain size He] is chosen 
as the critical length parameter. Thus, the local criterion is stated as; 


D = De at / = le. 


( 0 , 6 ) 


That is, whenever D reaches the critical value at a distance Ic ahead of the crack tip, 
the crack starts growing. 
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The crack growth initiation criterion cited above is used in numerical simulation 
to predict the fracture parameters like critical load (Per), fracture toughness [Jic] and 
critical crack-tip opening displacement (Stc) in plane strain condition. Standard fracture 
mechanics specimens ( CT and TPB ) are chosen for numerical simulation as well as 
for conducting experiments ( both mechanical and metallurgical ). The numerically 
predicted values of fracture parameters match favourably with their experimental values. 
Some additional results are also obtained from numerical simulations. These are also 
consistent with the results published earlier. Therefore, it is expected that the proposed 
local criterion ( equation (6) ) can predict the crack growth initiation in Mode-I ductile 
fracture. 

A three dimensional numerical simulation, using the above crack growth initiation 
criterion, is also attempted in the last part of the thesis. The results give a consistent 
qualitative trend. 

Finally, it is concluded that the critical damage Dc (or the critical area void fraction) 
is a satisfactory measure for predicting the crack growth initiation in ductile materials. 
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Chapter 1 


INTRODUCTION 


1.1 Introduction 

Application of linear elastic fracture mechanics to design of machine or structural com- 
ponents is restricted by the condition that the plastic zone at the crack tip be small 
compared to the size of the K dominant field and any relevant geometric dimension. It 
is virtually impossible to satisfy this condition for high toughness, low strength materials 
which generally undergo extensive plastic deformation and crack-tip blunting prior to the 
beginning of crack growth. Crack initiation in these materials is usually followed by a 
stable crack growth before it fails in catastrophic manner. The need to incorporate the 
influence of significant plastic deformation in crack growth initiation ( i.e. the onset of 
growth of a existing crack ) and subsequent stable growth is the main impetus for the 
development of elastic-plastic fracture mechanics. While the development in this field 
is ongoing, the theory of ductile fracture is not yet fully established. The use of HRR 
( Hutchinson, 1968 [1] , Rice and Rosengren, 1968 [2] ) stress field to predict the stress 
and displacement fields ahead of a crack tip is restricted by several limitations like pro- 
portional loading, small crack-tip geometry change etc. Therefore, the use of Jjc ( Rice, 
1968 [3] ) as a ductile fracture parameter becomes dependent on several conditions. 

Recent inv^tigations on ductile frcicture show that the crack growth initiation in a 
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ductile material is a microscopic phenomenon occurring over a small zone around a crack 
tip. This zone is termed as fracture process zone which is estimated to be of the order 
of crack tip opening displacement ((<5f)c) [ 4 ]. Within this fracture process zone, profuse 
void growth and coalescence occur which cause a local instability and ultimately the voids 
join in a void sheet to form micro cracks. When these micro cracks impinge with the 
parent crack, the crack growth initiation occurs. The void growth phenomenon depends 
on the purity of the material i.e. the presence of second phase particles or inclusions in 
the material. It also depends on the microstructural features of the material like grain 
size, inter particle spacing etc. For a material like steel the above phenomena are evident 
from the fractographs taken on the fracture surface where dimpled ruptured surface due 
to void growth is clearly visible. It is found from experiments that macroscopic material 
properties change when voids grow to a large extent. 

While formulating the problem of ductile fracture in the frame work of continuum 
mechanics, the major difficulty is how to represent the microstructural phenomenon 
with the help of continuum parameters. Here some of the most important issues are as 
follows. The first issue is about the constitutive modelling for incorporating the effects 
of void growth on the material behavior. The choice of appropriate micro-model for void 
coalescence and its representation in terms of continuum variables is another important 
aspect of the problem. Finally one has to make use of this coalescence model to formulate 
a criterion for crack growth initiation with the help of some local parameters. 

The above issues have been addressed by many research workers [ 5,6 ]. Due to 
difficulties in posing the problem in the frame work of continuum mechanics, several 
assumptions have been made. Thus, no single model addresses all the issues in a satis- 
factory manner. However, as far as constitutive modelling is concerned, the continuum 
damage mechanics (CDM) model [ 7,8 ] shows a lot of promise. The model has the capa- 
bility to take care of loss of load bearing capacity and strain softening resulting from void 
growth. Similarly, Thomason’s limit load model [ 9 ] seems to represent the phenomenon 
of void coalescence quite adequately. 
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Almost all the predictions made about ductile fracture have been qualitative. A 
quantitative prediction is rare. The present work aims at quantitative prediction of crack 
growth initiation in ductile materials. Elastic plastic finite element method coupled with 
CDM is employed for this purpose. Thomason’s model for void coalescence is used to 
propose a local criterion for crack growth initiation. Numerical predictions are compared 
with experimental results to validate the proposed criterion. 


1.2 Literature Survey 

1.2.1 Stress and Displacement Analysis at the Crack Tip 

Linear Elastic Fracture Mechanics ( LEFM ) 

In a linearly elastic material, for a nominally stationary crack subjected to tensile (Mode 
I ) opening, the local crack tip stresses and displacements can be characterised in terms 
oi K I singular field [ 10, 11 ]. Accordingly 

= -|L/«(«) + 0(rl) + ... (1.1) 

V27rr 

-=/,j(^) as r-^O. (1.2) 

y iTvr 

where, r is the distance ahead of crack tip, 9 is the polar angle measured from the crack 
plane, is a dimensionless function of 0, gi is a function of 9 and bulk modulus of the 
material and fi stands for shear modulus. The configuration independent amplitude factor 
Ki is called the Mode I stress intensity factor. This parameter uniquely characterises 
the local stress field ahead of a crack and can be used as an indicator for crack growth 
(Irwin,1957). This asymptotic stress field dominates the crack tip vicinity over a region 
which is large compared to the scale of micro structural deformation and fracture events 
involved. 
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One of the principal limitations of the above approach is that a state of small scale 
yielding must exist. From equation (1.2) it is apparent that as r tends to zero the stresses 
become unbounded at the crack tip. In reality, however, such stresses are limited by local 
crack tip yielding which occurs over a zone ahead of the crack tip known as plastic zone. 
A*' approximate estimate of the plastic zone for monotonically loaded crack can be taken 
as [ 12 ] 

for plane stress. (1-4) 

for plane strain (1.5) 

where ay is the yield strength of the material. Provided the extent of local plasticity is 
small compared to the extent of K / field which itself is small compared to the geometrical 
dimension of the body (including crack length) , the plastic zone can be considered as 
a small perturbation in the linear elastic field and the dominance of Kj field near crack 
tip remains valid. This situation is called the small scale yielding. 

Elastic-Plastic Fracture Mechanics 


1 1 


27r 


i.j 

f- 

Gtt 

Uy, 


Characterization of crack tip field by LEFM becomes unsuitable where small scale yield- 
ing conditions do not apply i.e when plastic zone size is comparable with the crack length. 
Since the Ki singular field is no longer appropriate in such instances, alternate asymp- 
totic analysis has been done to determine the crack tip stress and strain fields in the 
presence of more extensive local plasticity. Based on the deformation theory of plasticity, 
the asymptotic form of these local fields, for materials obeying the power hardening law 

( 1 , 6 ) 


is given by so called HRR [ 1, 2 ] singulax field: 

As r — )■ 0, 

J 

Otay^ylpf T' 


aij{r,9) 


(i.r) 
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Uj{r,6) aCyV 


J 


occrytylp^ r 

J 


N 

N+l 


^j{e,N), 


N 

WTT 


ujie.N) 


( 1 . 8 ) 


(1.9) 


lacTyeylN rj 

where a ajid N are the hardening parameters, cTe, and are the equivalent stress and 
plastic strain, CTy and Cy are the yield stress and yield strain, J;v is a numerical constant 
weakly dependent on N, and d'ij{&,N), eij(6,N), Uj{9,N) are normalised stress, strain 
and displacement functions of N depending on the condition of plane stress or plane 
strain. The amplitude of this field is J-integral ( Rice, 1968 [ 3 ] ). J uniquely charac- 
terises the elastic-plastic crack tip field provided some degree of hardening is present and 
thus can be used to predict crack extension. For small scale yielding J can be identified 
as the strain energy release rate 'G' and hence can be related to Kj through the relation 


J = G = ‘^ (1.10) 

where E' is the appropriate elastic modulus {E for plane stress, and Ej[l — v'^) for plane 
strain where v is the poisson’s ratio). 

It should be noted that HRR singular field and J-integral are strictly defined for 
a non lineax elastic solid for which stress is a function of current strain rather than for 
a physically more realistic elastic-incrementally-plastic solid where strain increment is a 
function of current state of stress. Provided the crack remains stationary and is subjected 
only to a monotonically increasing load, the plastic loading does not depart radically from 
proportionality and this field remains valid. However, for a growing crack where regions 
of elastic unloading and non proportional plastic loading are embedded within the J 
dominated field, the crack tip field is not properly modelled by the HRR field [ 13 ]. 

An alternative treatment of elastic-plastic crack growth initiation, which is not 
subjected to above restrictions, is to utilise the concept of crack tip opening displacement 
or CTOD [ 14 ]. From equation ( 1.9 ), it is apparent that the opening of crack faces 
varies as when r — »■ 0. This separation can be used to define CTOD (^j) as the 
opening where 45° lines emanating from the cradc tip intercept the crack faces [ 15 ]. 
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Thus, 


= dNi^y, N ) — for elastic plastic material. 


Kr 


ayE' 


for linearly elastic material. 


( 1 . 11 ) 


( 1 . 12 ) 


Here, d/v is a constant depending on N , Sy and the condition of plane stress or plane 
strain. For a perfectly plastic material and under plane stress condition, djv = 1. 


Effect of change of crack tip geometry 


Substantial progress has been made in the recent years in analysing crack tip stress and 
deformation fields in elastic plastic materials [ 16,17 ]. Most of the solutions are based on 
conventional ’small geometry change’ assumption. Two important changes are observed 
when the actual large geometry change of the crack tip is taken into account. Those 
are discussed by Rice and Johnson [ 18 ] . First, no severe strain concentration takes 
place directly ahead of the crack tip according to ‘small geometry change’ solution for 
nonhardening or moderately hardening materials. Further, an opening displacement at 
the crack tip can not defined. But, when progressive blunting of the crack tip is taken into 
account, intense strains do result directly ahead of the crack tip over a region of micro 
structural relevance like grain size, inter particle spacing etc. The second modification, 
due to large crack tip geometry change, has to do with stress distribution. Conventional 
’small geometry change’ solution gives a high stress concentration ahead of the crack 
tip ( infinite for a continuous hardening material as the crack tip is approached ) even 
when the strains are small. These high stresses result from high triaxiallity ( ^ ) at 
the crack tip. Here (Xm is the mean stress (cr^ = |<^fcit)- When crack tip blunting is 
taken into account and the proper boundary conditions are imposed at the crack tip, the 
triaxial.ity reduces substantially at the crack tip. In other words, the stresses at the crack 
tip become finite. The above observations are also found in finite element computation 
of McMeeking and Parks [ 16 ]. 

From the above discussion, it is clear that the intense strain at the crack tip due 
to crack tip blunting plays an important role in the ductile fracture process . Though 
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a strain based criteria is likely to be appropriate, still the role of triaxialaty can not be 
ignored. 


1.2.2 Classical Fracture Criteria 

The classical fracture criteria view fracture process in terms of continuum parameters. 
For brittle fracture, the crack tip stress field is characterised by Kj parameter whereas, 
for ductile fracture, it is J which specifies the crack tip stress and displacement fields. 
Sometimes crack tip opening displacement 5t is used to describe the crack tip field. Note 
that the crack tip opening displacement 8t has a unique relation with J or A'/ (equation 
1.11 and 1.12). Accordingly, the fracture toughness is defined by the critical value of 
K or J or 6f Thus, the crack growth initiation is predicted in Mode I by A'/ = Kic 
or J/ = J/c or 6t = ^tc- Normally, the critical value ( A'/^ or J/d is determined from 
experimental load deflection curve. ASTM ( E-399 ) specifies the standard for Kic and 
J/c tests. 

The critical value of fracture parameters can also be expressed in terms of a char- 
acteristic dimension of microstructure (Ic) using local failure models. In one such model, 
the critical crack tip opening displacement is expressed in terms of the mean distance 
between the void initiating particles (Rice and Johnson, 1970 [ 18 ]). Thus 

dtc = l to 2.7 X (1.13) 

where, Xo is the mean distance between the major void nucleating particles. From equa- 
tions ( 1.11 and 1.13 ), this model implies the fracture toughness as; 

J/c-^vX (1-14) 

although it is rare to find the toughness to increase directly with strength. Another 
approach to estimate the fracture toughness is to use the stress modified critical strain 
model ( Mackenzie and coworkers, 1977 [ 19 ] ). Here, local equivalent plastic strain has 
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to reach a critical ductility ( ), specific to the relevant triaxial ity { ^ ) , over 

some characteristic length Ic. Accordingly, the fracture toughness is found as: 

*J[c ^c* 

Unlike the critical CTOD criterion, this criterion shows that the fracture toughness is 
proportional to product of strength and ductility which is physically more realistic. 

1.2.3 Recent Development in Ductile Fracture 

Recent investigations of ductile fracture reveal that a ductile fracture occurs mainly due 
to microvoid nucleation, growth and finally coalescence into micro crack. It is found from 
experiments that macroscopic material properties change due to void growth. As a result, 
the constitutive equations need to be changed when the extent of void growth is large ( 
more than 10% ) [ 20 ]. Thus, a realistic model for prediction of ductile fracture should 
include a void growth dependent material behavior and a condition for void coalescence. 
Three broad approaches have emerged which try to predict ductile fracture { micro crack 
initiation ) on the basis of void nucleation. growth and coalescence. These are as follows: 

(i) Porous Plasticity model of Berg and Gurson [ 21, 22 ] 

(ii) Void Nucleation. Growth and Coalescence model (Goods and Brown- Rice and Tracy- 
Thomason Model) [ 28, 29, 26 ]. 

(iii) Continuum Damage Mechanics Model [ 7, 8 ]. 

These three approaches are discussed in the following paragraphs. Next, some phe- 
nomenological models of ductile fracture are discussed. When a micro crack initiation 
criterion is to be used for prediction of crack growth initiation, an additional ( local ) 
length parameter is needed. Some local criteria for crack growth initiation are discussed 
at the end of the section. 

Porous Plasticity Model of Berg and Gmrson 


In this model, the material with voids is idealised as a porous material. Thus, 
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its constitutive equation is derived from the plastic potential of porous material. In a 
porous plastic material, there is a possibility of dilatation because of which the yield 
surface does not remain an infinitely long cylinder but capped by elliptical surfaces. 
Based on Berg’s [ 21 ] model of dilational plasticity, Gurson [ 22 ] proposed a plastic 
potential for porous plastic material. Starting from a unit spheroidal cell with a single 
central void, he obtained the following expression for the potential function for a porous 
solid with randomly distributed voids of volume fraction ’/’ ; 

= T + / cosh _ 1 _ ,3/^=0 (1.16) 

Here 9i = 92 = 93 = 1- (Tvergaard [ 23 ] assinged some other values to the constants 
qi and 92 ( 9i = 1-5, ?2 = 1 and 93 = 9i ) to bring the predictions of the model into 
closer agreement with full numerical analysis of a periodic array of voids). S and are 
the effective and mean stresses of the porous aggregate and cTy is the flow-stress of the 
incompressible matrix. 

In this model, void nucleation is assumed to be a continuous process. Accordingly, 
it is represented by the relation 


nucleation — 

where A is constant. The void growth law is derived from the condition of plastic incom- 
pressibility. Thus 

= ( 1 . 18 ) 

The total incremental void volume fraction is of course the sum of the incremental void 
volume fraction due to nucleation and growth: 


df — df nucleation d" df growth' 

When the surrounding material ceases to deform, the dilation gets localised in a band 
(called as shear band). This concentration of dilational plastic deformation in a small 
zone gives rise to plastic instability leading to ductile fracture. Thus, in this model. 
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ductile fracture is regarded as the result of an instability in the weak dilatational plas- 
tic flow field in the shear band. Rudnicki and Rice [ 24 ] derived the conditions for 
the instability without considering any imperfection in the material. Later Tvergaard 
[ 23 ] and Yamamoto [ 25 ] investigated the conditions for instability in the presence 
of material imperfections. They represented the fracture criterion as a graph of critical 
localisation strain versus the critical void volume fraction with strain hardening exponent 
as a parameter. 

In this model, it is assumed that the void coalescence by internal necking of intervoid 
matrix is a secondary effect which can only develop after the incipient formation of an 
intense shear band. However, according to Thomason [ 26 ], the primary mechanism of 
ductile fracture is the internal necking of the intervoid matrix which can be adequately 
modelled only by the incorporation of a strong dilatational response. It is this failure to 
appreciate the dominant effect of internal necking, that severely limits the validity of the 
Berg-Gurson model. 

Before discussing the next model, it should be mentioned here that some research 
workers [ 27 ] have used the Berg-Gurson model without the instability conditions of 
Yamamoto. They characterised fracture by the critical value of ’/’ ( void volume fraction ) 
which was determined experimentally. 

Void Nucleation, Growth and Coalescence Model (Goods and Brown- Rice and Tracy- 
Thomason) 

In this model, Thomason [ 26 ] combined the results of Goods and Brown [ 28 ] 
on void nucleation, those of Rice and Tracy [ 29 ] on void growth and his own on void 
coalescence to arrive at a fracture criterion in the form of a graph of fracture strain 
versus the mean stress. For numerical predictions, this graph can be used along with 
conventional elastic plastic constitutive equations. Thus in the present model, the effect 
of void nucleation and growth are incorporated not in the constitutive equations but 
in the fracture criterion itself. The void nucleation model of Goods and Brown, void 
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growth model of Rice and Tracy and void coalescence model of Thomason are described 
below. Along with them, some additional models on void nucleation and gi'owth are also 
discussed. 

Void Nucleation Models 


A stress based void nucleation model is based on the condition that the void nucle- 
ation by decohesion of particles occurs whenever the normal stress at the particle/ matrix 
interface reaches the critical value ac- Goods and Brown [ 28 ] used this condition to 
derive the following relation for void nucleation strain (e") in terms of the mean stress 

e^,=K r{a,-cT„,? (1.20) 

Here, r is the particle radius and K is a material constant depending on the particle 
volume fraction. Experiments on Fe — Fe^C system [ 28 ] confirm the linear relationship 
between the void nucleating strain and the particle radius. The relation ( 1.20 ) is for 
small spheroidal particles (r < lum ). For larger size particles, the approximate analysis 
of Argon et al [ 30 ] gives the following condition for the void nucleation by decohesion: 

O-eg + <^m = (^c- ( 1 - 21 ) 

For a given problem, this condition can be expressed in terms of the nucleation strain. 
Note that this condition is independent of particle radius. 

A different void nucleation model has been proposed by Gurland [ 31 ]. His ex- 
periment on 1.05% C spherodised steel shows that the void nucleation by cracking of 
cementite particles takes place continuously. Voids nucleate at all strain levels depending 
on the size, shape and orientation with the maximum principal stress direction. It is 
observed that the volume fraction of broken particles is a linear function of the plastic 
strain. Note that the void nucleation relation (equation 1.17) used in the Berg-Gurson 
model is of this type only. 

It should be noted that the development of a very general model of void nucleation 
is a diflficult task as a typical commercial alloy contains a wide range of particle types 
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and particle morphologies which can result in a variety of void nucleation mechanisms 
operating simultaneously. Further such a model may be too complicated to be amenable 
to analysis for making useful numerical predictions. 

Void Growth Models 


Once the nucleation of a microvoid takes place either by decohesion or cracking of a 
second-phase particle or inclusion, the resulting stress-free surface of the void produces a 
localised stress and strain concentration in the surrounding plastically deforming matrix. 
With continuing plastic flow of the matrix, the void undergoes volumetric growth and 
change of shape. If it is assumed that the sites of void nucleation are sufficiently far 
apart to have any interaction between the neighbouring voids, then one can study the 
void growth phenomenon by considering a single void in an infinite medium. 


Rice and Tracy [ 29 ] considered a single spherical void of initial radius Ro in a 
remote uniform strain rate field e,j and remote stress field (7,j in a rigid-plastic material. 
They derived the following expression for the rate of change of the radii of cur\'ature [Rk] 
in the principal strain-rate directions: 


where 


Rk — [(1 -1- -t- teqD\ R„ 


(fc— 1,2,3) 


( 1 , 22 ) 


E 


Rme 


‘-eg 

1/ 


for linear hardening and for low values of am for non-hardening, 


D = 0.75 


1 for high values of am for non-hardening. 


for linear hardening, 


0.558 sinh f ) + O.OOSt' cosh 
V2 ayj 

^(-Rl + -^2 + Rs), 


( 3 (Tm 
\2 , 


e^q = equivalent strain rate, 


principal strain rates, 
3e2 


• • ? 
— £3 


Lode variable. 


for non-hardening, 
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Note that an initial spherical void grows into an ellipsoidal void of principal radii Ri , R 2 
and i? 3 . Thomason [ 26 ] integrated the above equation by assuming that the principal 
axes of the strain-rates remain fixed in direction throughout the strain path and obtained 
the following expressions for the principal radii of the void: 


where 


U 


B 


+ 3 , 


Bo, 







(1.24) 

(1.25) 

(1.26) 
(1.27) 


and ef is the integral of the largest principal strain rate. Thomason used the expressions 
(1.23 - 1.25) in his derivation of the condition for void coalescence. 


For an array of void nucleating particles of diameter Dp and spacing dp, setting the 
initial void radius as Dpjl and integrating equation ( 1.22 ) upto the point of fracture, 
the fracture strain tj ( for a non hardening material) is obtained as: 




‘“(fe) 


(1.28) 


0.28 exp ( 1-5^ ) 

Similar analysis has been done by McClintock [ 32 ] for the case of cylindrical voids in 
plane strain condition. The corresponding fracture strain expression is given by 


sinh (1 -n)(cra + CTi)/^ 


(1.29) 


Here CTa. and Cb are the principal stress components and n is the hardening parameter 
defined by 


^, = k{4,)\ 


(1.30) 
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The material constajat K is different than that used in equation (1.20). 

Some research workers [ 33,34 ] have used the expressions (1.28) and (1.29) for pre- 
diction of ductile fracture with the help of some experimentally determined parameters. 
The major limitations of this approach are that it ignores the effects of void nucleation 
( as the expressions are derived from a pre-existing finite size void ) and void coalescence. 
Thus, usually €/ is overestimated. However, these expressions do reveal the dependence 
of €/ and fracture toughness on triaxial.ity (^), hardening n and the purity of the 
material. 

Void Coalescence condition 

Thomason [ 26 ] has modelled void coalescence phenomenon as plastic instability 
due to necking of the intervoid matrix. According to him, the sufficient condition for 
plastic instability of the intervoid matrix is given by 

<TnA„-<7i=0 (1-31) 

where An is the area fraction of the intervoid matrix perpendicular to the direction of 
majcimum principal stress <ti and cr„ is the plastic constraint stress. He considered a 
geometrically equivalent square prismatic void with the same principal dimensions as the 
ellipsoidal void and used the upper bound method to obtain the expression for the plastic 
constraint stress in terms of the void dimensions, intervoid spacing and the yield stress. 
The details of this derivation are described in Chapter 2. 

Fracture Criterion 

Thomason used expressions (1.23 - 1.25 ) for the void dimensions to express the 
void coalescence condition ( equation 1.31 ) in terms of the void growth strain ef and the 
mean stress cr^n. By superposing this condition on the void nucleation relation ( equation 
1.20 ), he obtained the fracture criterion as a graph of fracture strain ( e/ = + ^ ) 

versus the mean stress. 

The main drawback of this approach is the use of expressions ( 1.23 - 1.25 ) for the 
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void growth. These expressions are based on integration procedure which assumes that 
the principal directions of strain rates remain fixed throughout the strain path. This is 
true only for the case of small strain and rotation. As a result, this approach can not be 
used when the strains and/or rotations are large. 

Continuum Damage Mechanics Models 

In continuum damage mechanics model, the change in material behavior due to 
void growth is taken care of by introducing the damage variable as an internal variable 
in the plastic potential. The damage variable is supposed to quantify the intensity of 
microvoids. As a result, it is usually identified as either the void volume fraction or area 
void fraction. The theory of continuum thermodynamics is used to derive the plastic flow 
rules and the damage growth law. 

It should be kept in mind that the main objective of continuum damage mechanics 
model is to incorporate the effects of void growth in the constitutive equations. However, 
there is a scope to include the void nucleation phenomenon in the damage growth law. 
But, as fax as void coalescence is concerned, it has to be incorporated as an additional 
condition in terms of the continuum parameters. This condition, which serves as a 
fracture criterion, has to be based on an appropriate micro model. Thomason’s limit 
load model [ 26 ] seems to be a good candidate for this purpose. 

Based on Kachanov’s idea of effective stress ( for creep damage ), Lemaitre [ 7 ] 
proposed a damage mechanics model for void growth in elastic plastic materials. Details 
of the model axe described in Chapter 2 of the thesis. Lemaitre used a simple damage 
growth law in which the damage rate is linearly dependent on the plastic strain and 
the thermodynamic force corresponding to the damage. Later, Tai and Yang [ 35 ] 
modified the damage growth law to make it dependent on the current value of damage. 
They obtained the material constants appearing in the damage growth law from the 
experimental results of LeRoy et al [ 36 ] on area void fraction measurements at differoat 
strain levels. Jun [ 37 ] introduced non linearity with respect to plastic strain in the 
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damage growth law. The degree of nonlinearity for a material is adjusted by a parameter 
which is determined from tension test data. Applications of Lemaitre’s model are found 
in the works of many authors, however, there seems to be no attempt at combining the 
model with a void coalescence criterion for the purpose of prediction of ductile fracture. 


Another popular model is due to Rousselier [ 8 ] . Rousselier model consists of a 
damaged plastic potential in terms of current density p of the material and mean stress 
CTm- Accordingly, the plastic potential is given as 


$ = 





(1.32) 


Here, is the deviatoric part of the stress tensor <7,j, J 2 is the 2nd invariant of cr[j, 
5 is a function related to the dialatational effect of void growth and B stands for the 
thermodynamic force corresponding to damage variable /? which is implicitly present in 
the plastic potential through the density function i.e 


» = m 


(1.33) 


The mass conservation law is used to derive the expressions for the function 5, the 
thermodynamic force B and the damage growth /?. 

Bilby, Howard and Li [ 38 ] applied the Rousselier’s model to predict the J — R 
curve for a spinning cylinder. The model was implemented with updated Lagrangian 
elastic plastic FE programme. However, for the purpose of prediction of crack growth, 
instead of using a condition for void coalescence, they used a stress based local criterion. 
Zheng, Luo and Zheng [ 39 | used Rousselier’s model to derive a macro damage parameter 
for prediction of crack initiation in metal forming processes. 

Phenomenological Models 

In the absence of reliable quantitative models for incorporating the phenomena of 
void nucleation, growth and coalescence in materials undergoing large plastic deforma- 
tion, in metal forming processes, empirical relations based on some phenomenological 
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models axe used. They are many in number. Here, only two popular models are dis- 
cussed. 


Cockcroft and Latham [ 40 ] suggested a fracture criterion of ductile material based 
on the idea of critical plastic work. Accordingly, their criterion is stated as 



aeqdtg^ = constant. 


(1.34) 


Here, Cf stands for fracture strain in tension test. For a tensile specimen, the change 
in neck geometry influences the fracture process. In order to take care of the effect of 
change in neck geometry, the above expression is modified as 



= constant 


(1.35) 


or, 



= constant 


(1.36) 


Here, cti is the maximum normal stress. The dimensionless factor acts as a stress 

concentrator. 


The main criticism of the model is that the role of hydrostatic stress is not incor- 
porated in this model even though it is accepted, in general, that the hydrostatic stress 
also influences the ductility of the material. 


Another populax criterion used in the field of metal forming is due to Oyane [ 41 ]. 
The Oyane’s model is based on a theory of plasticity of porous materials where dialata- 
tional stress-strain relation is assumed as 


de 


m — 


dCgq 

~A 



(1.37) 


where e^q is equivalent strain, stands for volumetric strain and A and Ao are material 
constants. Oyane integrated equation (1.37) over the total strain path to give 





g'm 

AoCeq 


dCeq. 


(1.38) 
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Here, Cj is the fracture strain in tension test. It is assumed that the left hand side of 
equation ( 1.38 ) is a material constant. Thus, 


i 


‘"i+ ■"- 


Act, 


dtgn — Bq, 


a constant. 


(1.39) 


eg, 


The main objection to this criterion is that the ductile fracture is governed by criti- 
cal volumetric strain independently of deviatoric strain. Ductile fracture in pure shear 
contradicts this theory. 


Local Crack Growth Initiation Criteria 

All the criteria discussed so far are for micro-crack initiation. When these criteria 
are to be used for prediction of crack growth initiation, an additional ( local ) length 
parameter is needed. Some existing local criteria for crack growth initiation are discussed 
next. 


Critical Stress Criterion 

A stress based criterion states that the cleavage fracture propagates in an unstable 
manner when the maximum principal stress {(Tyy) ahead of a crack tip reaches a critical 
value RKR ( Ritchie, Knott and Rice ) [ 42 ] proposed a local criteria for cleavage 
fracture of a sh«irp crack which requires the local tensile stress {cTyy) to exceed a critical 
fracture stress over a microstructural characteristic length Ic- In their work on quantita- 
tive prediction of fracture toughness of mild steel at low temperature, the fracture stress 
was determined by breaking a V notched bent bar at low temperature and a-yy was found 
from Hill’s slip line solution. RKR found that their prediction gives best agreement with 
the experimental value when Ig = 2<i, where d is the diameter of Fe grain. 

Critical Strain Criterion 

Ritchie, Server and Wullaert [ 43 ] used the strain based criterion to predict fracture 
toughness of A533B and A302B alloy steels. The critical fracture strain was determined 
from tension test of circumferentially notched round specimen. The equivalent plastic 
strain ( ef,) was found from Rice and Johnson’s [ 18 ] solution. Their findings show a 
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critical length Ic, varying from one to six times the average MnS particle spacing. 


1.3 Objective of the Present Work 

The present work is an investigation of Mode I fracture in ductile material where void 
nucleation, growth and coalescence play a predominant role. The aim of the work is to 
identify a critical continuum parameter and critical (local) length parameter to predict 
crack growth initiation in Mode I ductile fracture. The following are the broad objectives 
of the thesis. 

1. To study the role of void growth in ductile fracture process and to propose a 
void growth law. 

2. To determine the critical value of area void fraction of the material for micro 
crack initiation and to establish its independence of specimen geometry. 

3. To propose a crack growth initiation criterion in Mode I ductile fracture using 
the critical area void fraction and a critical ( local) length parameter. 

4. To predict the global fracture paxameters (like .7/^ ) numerically using the 
proposed crack growth initiation criterion 

5. To verify the numerically predicted values of the fracture parameters with ex- 
periments ( both mechanical and metallurgical ) to validate the proposed crack growth 
initiation criterion. 


1.4 Plan of the Thesis 

Lemaitre’s continuum dannage mechanics model is used to derive the constitutive equa- 
tions (flow rules and damage growth law ■) of a material with voids. The damage is 
identified as an area void fraction. The updated Lagrangian elastic plastic FE formula- 
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tion ( which includes both material and geometric nonlinearities ) is used for numerical 
simulation. Thomason’s criterion for void coalescence is used to determine the critical 
value of the axea void fraction. The critical value of area void fraction along with a 
critical length parameter is used to predict the global fracture parameters. Experiments 
have been carried out to verify the results of numerical simulation. 

Chapter 2 of the thesis deals with constitutive equations of a material with voids. 
The concept of continuum damage, due to Lemaitre[ 7 ], is used to modify the von Mises 
plastic potential. A damage growth law, based on LeRoy’s experimental results [ 36 ] 
on AISI 1090 spherodised steel, is proposed. The law also incorporates void nucleation. 
A micro model for void coalescence ( Thomason [ 26 ] ) is also discussed. The model is 
modified appropriately, without changing the basic principle, to suit the elastic plastic 
FE formulation. 

Chapter 3 is devoted to updated Lagrangian elastic plastic FE formulation. Both 
the nonlinearities, material and geometrical, have been incorporated in the FE analysis. 
Deformation due to finite strain as well as rotation has been considered. Different nu- 
merical schemes are also discussed here. The elastic plastic FE programme developed on 
the basis of the above formulation is used in the later part of the thesis to investigate 
the Mode I fracture of ductile material numerically. 

Chapter 4 is a numerical investigation of micro crack initiation, due to void nucle- 
ation, growth and coalescence, for different specimens of AISI 1090 spherodised steel. 
The constitutive equations developed in chapter 2 and the FE formulation of chapter 3 
axe applied together along with Thomason’s model. The critical value of damage variable 
(area void fraction) is found to be almost constant for all the specimens analysed. Hence, 
the critical value of damage variable is taken as a paxameter for crack growth initiation 
in the later part of the thesis. 

Chapter 5 deals with the Mode I fracture of a ductile material like AISI 1095 
spherodised steel in plane strain condition. A crack growth initiation criterion is proposed 
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in terms of the damage as a critical continuum parameter and austenite grain size as a 
critical length parameter. Numerical simulations of standard fracture mechanics tests 
have been carried out using the FE formulation of chapter 3. The results of numerical 
simulations are compared with experiments ( both mechanical and metallurgical ) to 
validate the proposed crack growth initiation criterion. 

Chapter 6 is an extension of the work to 3D case. The effects of variation in thickness 
on critical fracture parameters are studied. The variation of continuum parameters at 
the crack front, along the thickness, has also been observed in this study. 



Chapter 2 


CONTINUUM DAMAGE 
MECHANICS 


2.1 Introduction 

Investigation of ductile failure of metals requires the knowledge of material behavior at 
microstructural level. It has been found from metallurgical test results that the ductile 
fracture occurs due to void nucleation, growth and finally coalescence into a micro crack. 
In order to model the void growth phenomenon in continuum mechanics, the idea of 
continuum damage mechanics (CDM) is very helpful. In CDM, the material behavior is 
represented by a plastic potential which includes damage as an internal variable. The 
plastic flow-rules and the damage growth law are derived from the plastic potential. 

The basis of the continuum damage mechanics rests on the theory of continuum 
thermodynamics. It is well known that the constitutive relations must obey the basic 
laws of thermodynamics. Hence a brief discussion of continuum thermodynamics is given 
first in sections 2.2 and 2.3 [ 44,45 ]. Next, Lemaitre’s continuum damage mechanics 
model and the resulting constitutive equations of a damaged material are presented in 
sections 2.4 and 2.5. In the end, a criterion for void coalescence and micro crack initiation 
is discussed. 
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2.2 Thermodynamics of a Continuum 


A thermodynamic process is normally described by a set of kinematic variables cik {k = 
1,2, ... ,n) and the absolute temperature 9 (non negative). These quantities are known 
as independent state variables. Any function of the state variables is referred as state 
function. Examples of state function are internal energy U{ai:,9), entropy S{ak-9), free 
energy 0) etc. For a reversible process, the free energy or the thermodynamic 

potential completely specifies the thermomechanical behavior of a continuum. By means 
of the two thermodynamic laws, we get the following expression for the conservative part, 
Al, of the thermodynamic force corresponding to a^: 


Ai 


dak' 


( 2 . 1 ) 


Ak is also called as the conjugate variable corresponding to a^.. Since, the entropy pro- 
duction is zero, there is no need to introduce the dissipative part of the thermodynamic 
forces. As an example, consider the following expression for specific free energy for an 
elastic process: 


i> = Y 


Os. 


( 2 . 2 ) 


Here, the components of elastic strain tensor are the kinematic variables. Cijki is the 
fourth order elasticity tensor and s is the specific entropy. The Cauchy stress tensor is 
given by 


(Tij = p 


d'lj) 


Cijki^ki' 


(2.3) 


For an irreversible process, it is necessary to consider the dissipative thermodynamic 
forces in order to describe the thermomechanical behavior of a continuum. For such 
a process, the rate of entropy production is not a state variable but in general a 
function of djt, akj and 9. Normally this functional relationship is expressed for the 
dissipative power which is nothing but the product of absolute temperature $ and the 
rate of entropy production 5*. The superscript ’ i ’ stands for the irreversible process 
while the superscript denotes the time rate of the quantity. Thus for an irreversible 
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process, the dissipative power can be expressed as 

6S' = $( fljk, 6, ak^ 9 ). 


(2.4) 


If we assume that an irreversible system is purely dissipative, then the dissipative forces 
Af can be completely derived from Thus 




(2.5) 


dak' 

The quantity $ is called the dissipative potential. The 2nd law of thermodynamics states 


that 


$ > 0 . 


(2.6) 


2.3 Thermodynamic Process with Internal Variables 


In the field of mechanics or thermodynamics, all variables can be classified into two 
categories, namely the variables which are measurable or conti'ollable and the variables 
which can not be measured or controlled directly. The first kind is termed as external 
vaxiable whereas the second one is called as internal variable. When we include them 
in thermodynamic potential, we call them as external state variable and internal state 
variable respectively. 


For dissipative phenomena, the current state also depends on the past history. This 
history cm only be represented by means of internal variables. For example in a damaged 
elastic-plastic material, the history of deformation can be represented by accumulated 
plastic strain, accumulated damage etc. 


The thermodynamic potential or free energy as a function of internal state variables 
Ojk is expressed as 


^ = ^(afc,a]t,^) 


(2.7) 


and the corresponding thermodynamics forces are given by 


Ai = 


dak 


( 2 . 8 ) 
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r- 

dak' 

(2.9) 


(2.10) 

do' 


Here, is the conservative part of the thermodynamics force (3^ corresponding to the 
internal state variable at. 


When the dissipative potential $(dfc, o:*;, at, 0 :^, 9) is also function of internal state 
variables and dfc, the dissipative forces are given by 

dcLk 


Ai = 


( 2 . 11 ) 


and, 



( 2 . 12 ) 


Here, is the dissipative part of the thermodynamic force /?fc corresponding to the 
internal variable ak. Since /5jt is an internal force. 


ffk + Pt^O 

or, 

= -01- (2.13) 


2.3.1 Elastic-plastic Process 

When the temperature change is not significant, the dissipative potential for an elastic- 
plastic process is expressed as 

* = «(er^,p,i),e?,p,n). (2.14) 

Here, p and D are taken as internal variables corresponding to isotropic hardening and 
damage respectively. The physical identification of the damage variable is discussed in 
detail in the next section. For the case of strain hardening, p is nothing but the equivalent 
plastic strain Thus 

P = <, = / 

Jo 


(2.15) 
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and 



(2.16) 


The quantity is the plastic part of the strain rate tensor and efj is the integral of efj 
along the particle path. Therefore, 


o-.i = 


K' 


(2.17) 




(2.18) 

(2.19) 


Here, —R and —Y are the dissipative parts of the thermodynamic forces corresponding 
to the internal variables p and D respectively. These laws are called as complementary 
laws while equations ( 2.8 - 2.10 ) are called as state laws. Using Legendre-Fenchel 
transformation [45], one can transform p, !)■, efj, p, D) to its dual —R-, —Y) 

and write the complementary laws as the evolution laws of the flux variables. Thus 




da-ij ’ 




p = 

D = 


d{-Ry 

di-YY 


( 2 . 20 ) 

{ 2 . 21 ) 

( 2 . 22 ) 


Since or $ is governed by an inequality ( see equation 2.6 ), its actual value 
is usually difficult to determine. As a result, the above laws are not useful from the 
computational point of view. To derive the useful relations, we maximise the dissipation 
using the yield function or the plastic potential F(crij,—R,—Y), which has a zero 
value on the yield surface, as a constraint. Thus 


-\F) = 0 


(2.23) 


where 6 denotes the variation and A aots as a Lagrangian multiplier. The function F is 
called an indicator function. Combining equations (2.20 - 2.22) with equation (2.23), we 



get 
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dcTii’ 

(2.24) 

OF 

di-RY 

(2.25) 

f)F 

a(-Yy 

(2.26) 


The laws (2.24 - 2.26 ) are called as plastic flow rules and they describe the mechanical 
behavior of an elastic-plastic material. 


2.4 Introduction to Damage Mechanics 


Damage represents surface discontinuities in the form of micro cracks or volume discon- 
tinuities in the form of micro voids. The description of material behavior of a damaged 
material involves an additional internal variable, called as the damage variable, which 
quantifies the intensity of damage. If it is assumed that the voids are scattered in an 
isotropic way, then this variable can be assumed as a scalar quantity which is denoted 
by D in the literature. Here the damage variable is identified as area void fraction at a 
point in a plane. That is. 


D = 


AA„ 

AA 


(2.27) 


where AA is an infinitesimal area around the point in some plane and AA,, is the area 
of the void traces in the plane contained in AA. The introduction of damage variable 
leads to the concept of effective stress i.e. the stress calculated over the effective area 
( A A — AA„ ) that actually resists the forces. Thus the effective Cauchy stress tensor 
at a point is defined as 


rr* — 


(2.28) 


Another important concept is the principle of strain equivalence which states that the 
deformation behavior of a damaged material is represented by the constitutive laws of 
the virgin material in which the usual stress is replaced by the effective stress . The 
concept is represented in Fig( 2.1 ). 
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When temperature change is not significant, the specific free energy of an elastic 
process in a damaged material can be obtained using the principle of strain equivalence. 
Following Lemaitre [ 7 ], we get 

V- = YCiik,e%<iUl - D). (2.29) 

From equation ( 2.29 ) the Cauchy stress components are found as 

=l>^ = Ciik,4,(l - D). (2.30) 

Further, the conservative part of the thermodynamic force (F) corresponding to D is 
given by 

r = pH = (2-3i) 

Since it is an internal force, the dissipative counterpart of Y is obtained as ( see equation 
2.13 ) 

= (2.32) 

For a better physical interpretation of Y, equation (2.32) can be modified as follows. 
From equation (2.30), we get 



dD ~ dD 

(2.33) 

For constant stress. 

II 

O 

(2.34) 

and hence. 

C<iu4, = C„«(l - D)^. 

(2.35) 


Using equation ( 2.35 ) in equation ( 2.32 ), we get 


-Y = 


- D)4i^ 


(T^constant 


or, 

_ 1 _ 1 dW, 

2 dD 


a^constani 


(2.36) 
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Here dW^ = (Tijde^j is the elemental strain energy per unit volume. Therefore, —Y is 
nothing but the change in strain energy due to damage when the stress remains constant. 

An explicit expression for (— i^), for an isotropic material, can be found from equa- 
tions (2.32) and (2.30) which will be useful in the later part of the thesis. First, inverting 
the equation (2.30), we get 


1 


rSijkicTki 


-tj ~ j (2.37) 

where Sijki is the inverse of Cijkl- The above expression, when substituted in equation 
(2.32), gives —Y in terms of the stress tensor. 

1 


Y = 


;Sijkicrijcrki- 


(2.38) 


2{1-Dy 

Using Sijki for an isotropic material and decomposing stress tensor in its deviatoric and 
mean components, we can write 

1 


Y 




E (1-T»)2 


E ' 


Here, E and i/ axe the elastic constants of the material. 


is the deviatoric part of cTij and 


, 1 
O’tj — Cfij — 


O’m — rCTii 


is the mean part. Using the definition of equivalent stress 


O’eg — 


3 , , 
L2^«i^u 


the expression ( 2.39 ) can be written as 


-Y 


-eq 


|(l + </) + 3(l-2.-)(^)^ 


'eq J 


(2.39) 


(2.40) 


(2.41) 


(2.42) 


(2.43) 


2E{l~DY [3' 

The term within the bracket in equation ( 2.43 ) is termed as triaxial’Ity function / (~^) 
whereas the quantity (^) is named as triaxialiity. Thus, 


-Y = 


■eq 


f 


2E{1-Dy-' \cr, 


(2.44) 


69. 
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Using the definition of effective stress, we get 

^*2 / 

-Y=^f h 

2E-' ' - 


(2.45) 


eq , 


2.5 The Constitutive Equations for Elastic-plastic 
Process in a Damaged Material 


The constitutive equations for elastic-plastic behavior of a damaged material can be 
derived from the appropriate plastic potential F. We assume that, in metal plasticity, it 
is possible to decompose F as 

F = + Fd{-Y,p,D) (2.46) 

where Fd is the plastic potential associated with damage such that it reduces to zero 
whenever D = 0. For a material yielding according to von Mises criterion, the form of 
Fi is 

— R 

(To (2.47) 




1-D 

where <To is the initial yield stress in tension. When D = 0, Fi reduces to the original 
form due to von Mises. Now the plastic flow rules ( equation 2.24 - 2.26 ) become 




dFi 

da; 


A 3 cr-j 


V 


eg 


p = A 


1 - D 2 <T, 

dFi _ A 
d{-R) ~ 1-D' 
OFd 


^ h{-Yy 

Combining equations (2.15) and (2.49), we get the following expression for A: 

A = (1 - D)eL. 


(2.48) 

(2.49) 

(2.50) 

(2.51) 


Substitution of expression (2.51) into equations (2.48) and (2.50) leads to the following 
constitutive equations ( stress-strain rate relation and damage growth law ): 


(2.52) 



31 


= (1 - (2-53) 

This stress-strain rate relation is not really convenient for an updated Lagrangian formu- 
lation for which the incremental stress-strain relation is needed. It is derived in Chapter 
3 starting from the expressions (2.46) and (2.47) for the plastic potential. 

Unlike i^i, Fd is not well established in literature. As a result, we can not use 
equation (2.53) as the damage growth law. Therefore, experimental results on void 
measurement at different deformation levels are used to propose a damage growth law. 
Based on experimental results of LeRoy et al [ 36 ] for spherodized steel, the following 
growth law is proposed. 

tl = ci%+(a, + a^D){-Y)il,. (2.54) 

Here the coefficients ai.aj and c are material constants which axe evaluated from ex- 
perimental results. The 1st term of the equation (2.54) represents the void nucleation 
as proposed by Gurland [ 31 ]. As stated in chapter 1, this model assumes that void 
nucleation takes place continuously at all strain levels. The other terms of the equation 
represent the void growth phenomena as proposed by several workers [ 35,37 ]. It is wor- 
thy to note that in Lemaitre’s [ 7 ] analysis the term corresponding to Oj is not considered 
where as in the work of Tai and Young [ 35 ] the term corresponding to ci is missing. 
From the experimental results of LeRoy et al [ 36 ] , it is observed that the graph of area 
void fraction versus strain is almost linear at low strain level but it is highly nonlinear 
at higher values of strain. Therefore it seems to be wise to retain both the terms in 
void growth rate equation ( 2.54 ). One can add more nonlinear terms like asD^, 04!?®, 

etc. But it is found that, for the material under study, such terms have insignificant 

\ 

contribution. 

Coming back to equation ( 2.45 ), the conjugate vaxiable(— F) corresponding to 
damage D can be expressed in terms of equivalent plastic strain in the following way. 
Using the definition of effective stress, the expression (2.47) for the plastic potential Fi 
can be expressed as 

Fl = < - (Ty 


(2.55) 
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where the yield stress ay of the material is given by 

~ 2 _ £) "*■ ^0^ • (2.56) 

For ' lost tt-' -^als, the c' r- b .'-dining relationship is v..- - by a power law: 

IT, = K (2.57) 

where K and n are hardening parameter and exponent respectively. Since Fi is zero 
along the yield surface, from equation (2.55) and (2.57). we get 

<, = !< {<,)"■ ( 2 -®) 


After substituting this expression in equation (2.45), —Y becomes 



(2.59) 


Then the damage growth law ( equation 2.54 ), in terms of the equivalent plastic strain, 
can be written as 


b = < + (ai + a,D)^ f il,. (2.60) 

The form of the damage growth law as cited in equation( 2.60 ) is used in FE computation. 


2.6 The Void Coalescence and Micro- crack Initia- 
tion 

Experimental studies on micro structural features of ductile fracture indicate that the 
ductile fracture process is essentially a localized plastic instability occurring simultane- 
ously in the inter void matrix between a large number of coalescing microvoids [ 26 ]. 
According to this reference, the sufficient condition for plastic instability of the intervoid 
matrix at a point is given by ( equation 1.31 ) 


ai = a„A„ 


( 2 . 61 ) 
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where <Ti is the current maximum principal stress of the macroscopically homogeneous 
state of stress at the point, An is the current area fraction of the inter void matrix 
perpendicular to the direction of ax and an is the critical stress required to initiate 
localized plastic flow or internal necking. The stress an is called the plastic constraint 
stress. 


Since there are about 10® microvoids in an area of 0.25 mm? of the ductile fracture 
surface[ 26 ], one can describe the fracture process by using a simple unit cell model in 
which the unit cell represents the statistical average of the microvoid size and intervoid 
spacing at the point under consideration. Thomason [ 26 ] considered a geometrically 
equivalent square prismatic void with the same principal dimensions as the ellipsoidal 
void ( Figs 2.2(a), (b) ) and used the upper bound method to obtain the expression for 
plastic constraint stress. A brief outline of the method is given in the next few pax’agraphs. 

Thomason assumed a certain portion of the unit cell as a plastic zone and divided 
that zone into a number of generating segments. He used parallel generating segment 
( Fig 2.2 (c) ) for a smaller void size and triangular generating segment ( Fig 2.2 (d) ) 
for larger void size. He assumed the following kinematically admissible velocity field. 


(a) For parallel generating segment ( Fig 2.2 (c) ): 

W 


u = 


V — 


2ax 

Wy 

2ax^ 


[{& + dY — , 

[{b + df - , 

Wz 


w = 


a 


s = \/u^ + 

(b) For triangular generating segments ( Fig 2.2 (d) ): 


Wd ,,, „ , 

U = — [(6 + d) + x] , 


V 


2ax 

Wyd 

2ax^ 


[(6 + d) + x] , 


(2.62) 

(2.63) 

(2.64) 

(2.65) 

( 2 . 66 ) 

(2.67) 
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Wzd 
‘lax ’ 


^ = ^[ib + d) + x] 




( 2 . 68 ) 

(2.69) 


Here u, v, w axe Cartesian velocity components at a point (x, y, r), is the velocity ( in 
z direction ) of the intervoid matrix at the top and bottom boundaries, a and b are the 
void dimensions, d is the intervoid spacing and i is the tangential velocity discontinuity 
on boundaries between the plastic zone and the remaining rigid zone ( such surfaces are 
indicated in the Figs 2.2 (c) and 2.2 (d) ). 


The rate of internal energy dissipation 7 in the plastic zone V' and the rigid-plastic 
boundaries S is given by 







(2.70) 


Here, a* stands for the yield stress of the rigid perfectly plastic material in the intervoid 
matrix. Using the strain rate velocity relations 


dii 


'-yy 


dv 


dw 

w 


\ (dii dv\ 

_ 1 / dw\ 

2 j ’ 

1 / dw du\ 

2 (&"'■&) ’ 


(2.71) 


and the expressions ( 2.62 - 2.69 ) for the velocity field, I can be expressed in terms of 
W, the void dimensions ( a and b ) and the intervoid spacing d. 


The rate of external work done is given by 


E = 


(2.72) 
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Here An represents the actual cross sectional area of the inter void matrix perpendicular 
to W. Setting £ = /, aji expression for an can be obtained in terms of the current void 
dimensions and <t*. Thomason obtained the solution of the above equation numerically. 
The numerical result obtained by him can be approximated by the following relation: 



(2.73) 


For a strain hardening material, the above expression should be modified by replacing 
a* with the current yield stress ay of the intervoid matrix material. Thus. 



(2.74) 


Elimination of cr„ from equations ( 2.61 ) and ( 2.74 ) leads to the following criterion 
for the void coalescence and micro crack initiation. 


ai - 



(2.75) 


Therefore, whenever the combination of stress components and void size and spacing 
satisfies equation ( 2.75 ), void coalescence and micro crack initiation take place at that 
point. Since the upper bound solution gives an overestimated plastic constraint stress 
(cr„), the above combination is an overestimation of onset of micro crack initiation. 


In order to incorporate the above criteria in FE analysis, the current void dimensions 
and spacing should be expressed in terms of stress or strain. Thomason used Rice and 
Tracey’s[ 29 ] expressions (equations 1.23 - 1.25 ) to express a, b and d in terms of strain 
and triaxiallity. However, as stated in Chapter 1, these expressions are valid only for the 
case of small straiin/rotation and therefore not suitable for large strain/rotation. Several 
research workers [ 46 ] have suggested that, at void coalescence, the a jd ratio should lie 
between 0.8 to 1.2 depending on the triaxial ity at the point. In the experimental 
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results of LeRoy et al [ 36 ], a/d is observed to be close to 1.0. Therefore, in the present 
work , it is taken as 1.0. Thus, 

^ = 1.0. (2.76) 

To relate the void dimension ' h ' and the intervoid spacing ' d' to the strain at the point, 
we assume that the equivalent plastic strain at the point, as defined by equations 
( 2.15 - 2.16 ), is more or less equal to the axial plastic strain of the unit cell ( Fig 2.3 ). 
Then, neglecting the elastic strain, we get 



An = ^ = exp{-el^) 

(2.77) 

and hence, 



(2.78) 

Equation ( 2.78 ) imphes that 



(2.79) 

Substituting equations ( 2.76 ^ 

), ( 2.77 ) and ( 2.79 ) into equation ( 2.75 ' 

1, we get 



Thus, whenever equation ( 2.80 ) is satisfied in the course of FE incremental computation, 
the micro crack initiation or void coalescence takes place at that point. 




V <7 \ a \ cr* 

Virgin Material Damaged Material Eqni\nlent virgin 

Material 


Fig. 2.1 : Strain equivalence principle. 





Fig. 2.2 The unit cell and velocity fields in 
the intervoid matrix. 

(a) An ellipsoidal void, (b) The equivalent square-prismatic 
void, (c) The parallel velocity field, (d) The triangular ve- 
locity field ( after P. F. Thomson [26] ). 









Chapter 3 

FINITE ELEMENT 
FORMULATION 


3.1 Introduction 

In this chapter a nonlinear finite element formulation has been developed using updated 
Lagrangian approach. An elastic-plastic matrix relating incremental stress and strain 
is derived using the plcistic potential proposed in Chapter 2 for a damaged material. 
Geometric non linear relationship between strain and displacement is introduced in the 
finite element formulation following a procedure given by Bathe [ 47 ]. Finally the 
numerical schemes, required to implement the finite element formulation, are discussed. 


3.2 Introduction to Non-linear Analysis 

In linear elastic finite element formulation, we assume that the displacements and their 
derivatives are infinitesimally small and the material is linearly elastic. The former as- 
sumption leaxis to the condition that the strain displacement matrix {B\ of each element 
is independent of nodal displacements. Moreover, the stress strain matrix [C] also re- 
mains constant for a linearly elastic material over the entire loading. The finite element 
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equilibrium equations become 

[I<]IL=F ( 3 . 1 ) 

where [A ] is the global stiffness matrix and U_ and F_ are the global displacement and 
load vectors respectively. The global stiffness matrix is an assemblage of the element 
stiffness matrices, i.e. 

[K] = UK.] ( 3 . 2 ) 

e 

where 

[IQ = J^[Bf[C][B]dv. (3.3) 

Here, V is the volume of the element and the summation is to be carried out after 
expanding the matrices [A^e]* Since the matrix equation (3.1) is a set of linear equations, 
the displacements are linear functions of load. 

In nonlinear finite element analysis, we encounter two types of nonlinearities , 
namely the material nonlinearity and the geometric nonlinearity. The 1st kind deals 
with the nonlinearity in material stress strain relation. Thus \C\ matrix becomes a 
function of the current state of stress. The other type of nonlinearity comes from the 
nonlinear relation between the strain and the displacement which is predominant in large 
deformation and rotation. Here the strain displacement matrix [5] becomes dependent 
on displacement. Thus different problems can be categorised in the following manner. 

Nature of the problem 
{a)Elastic analysis with small strain 

{b)Elastic analysis with large deformation 
and rotation. 

[c) Elastic — plastic with small strain Materially nonlinear only. 

{d)Elastic — plastic with large deformation Nonlinear in both material 
and rotation and geometric relations. 


Catagory 

Linear in both material 
and geometric relation. 

Geometrically nonlinear only. 
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3.3 Elastic-Plastic Analysis with Large Deforma- 
tion /Rotation 

3.3.1 Incremental Stress-strain Relationship representing Ma- 
terial Nonlinearity 


Since the stress-strain relationship is path dependent for an elastic-plastic material, it is 
convenient to represent it in incremental form. The incremental stress-strain relation is 
expressed as 

da=[C^^\de (3.4) 

where dgi and de are the incremental stress and strain written in vector form. The 
expression for the elastic-plastic matrix is derived from the plastic potential F. 

Initially, we shall derive the expression for for a material without damage. Since 

D is zero for such materials, the expressions (2.46) and (2.47) for the plastic potential 
reduce to 

F = Fi = (Teg — CTy (3.5) 


where now the yield stress is given by 


(Ty = {R + (To). (3-6) 

Note that F depends on the Cauchy stress tensor cr.j through cr^q and the isotropic 
hardening variable p through (Xy ( the variable p has been defined in the previous chapter 
by equations 2.15 - 2.16 ). For the derivation of [C®^] matrix, it is convenient to express 


(Xij in a vector form a: 


= {<Xu,CT22,...}- 


(3.7) 


Now the (von Mises) yield condition can be expressed as 

F( 21, p ) = o-eg( 21 ) - P ) = 0- 


(3.8) 


On the yield surface, dF = 0 and hence 

©' 


da+ -^dp = 0, 
op 


(3.9) 
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or 


crda — AXdt = 0 


(3.10) 


where, the flow vector a is defined by 


a = 


da' 


the parameter A is given by 


A = — T-r^i^dp. 


(3.11) 


(3.12) 


Xdt dp 

and A is the same scalar which appears in equations ( 2.24 - 2.26 ). When D = 0 { i.e. 
when F = Fi), equations (2.48) and (2.51) give the following expressions for the flow 
vector a and the scalar A : 

(3.1.3) 


" 2<7e,-’ 


Here a is the vector form of the deviatoric part of aij. Equation (2.15) gives 




(3.14) 


(3.15) 


Substitution of equations (3.8), (3.14) and (3.15) in (3.12) leads to the following expres- 
sion for A 

rJrr Arr 

( 3 . 16 ) 


^ 


dp delq 

The total strain increment can be split into elastic and plastic parts. Thus, 


de — d§^ + dF 
= [C]-^da + Xdt 


da 

= [C]~^da + Xdt a. 


(3.17) 


Here, [C\ is the matrix form of the elasticity tensor Cijki ( see equation 2.3 ). The 2nd 
paxt of the right hand side of equation( 3.17 ) is due to the flow rule (equation 2.24 ). 
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Pre-multiplying both sides of equation ( 3.17 ) by a^[C] and using equation ( 3.10 ) one 


can get 


Thus, 


which leads to 


\dt 


[C] de 

A + [C] a 


dt=[C]-^da-\- 


^ [C] de 
A + [C] a 


(3.18) 


(3.19) 


Therefore, 


d(j = 


[C]- 


[C^^] = 


[C] 


[C] a ^ [C] 

A 4- [C] a\ 

[C] a ^ [C] 


de. 


(3.20) 


( 3 . 21 ) 


A + [C] a J 

where the definitions of [C], a and A are as before. The expression for a is given by 
equation (3.13). For a material with power strain hardening law ( equation 'l.ol ), the 
expression (3.16) for A becomes 


A=.Kn (e;,) 


n— 1 


( 3 . 22 ) 


For an isotropic material, the expression for [C] for the 3-D 

case 

is given 

by 


‘1-1/ 

V 

V 

0 

0 

0 ■ 



V 

l-v 

u 

0 

0 

0 


j. , E 

u 

V 

l-u 

0 

0 

0 


(l + i/)(l-2i/) 

0 

0 

0 

l-2i/ 

2 

0 

0 



0 

0 

0 

0 

l-2t/ 

2 

0 



0 

0 

0 

0 

0 

l-2t/ 

2 - 



(3.23) 


From this, one can easily extract the expression for either the axisymmetric case or the 
plane strain case. 


The tensor form of the incremental relation (3.4) can be expressed as 


At,,. = C?£(i£« 


(3.24) 
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where da{j is the incremental stress tensor, tki is the incremental infinitesimal strain 
tensor and is the ( fourth order ) tensor form of [C^^\ matrix. The tensor dcr,j 
should be an objective stress tensor i.e. it should not be affected by a rigid body motion. 
Therefore daij can not be interpreted as an incremental Cauchy stress tensor. Two 
commonly used objective stress tensor in the literature are the 2nd Piola-Kirchoff stress 
tensor S'.y and the tensor based on Jaumann stress rate aij. The definitions of the tensors 
are given in the following sections. Thus daij should either be dSij or aijdt ( Here dt is 
the incremental time ). 


For a damaged material, the plastic potential is given by equations ( 2.46 - 2.47 ). 
When a similar derivation is carried out starting from equations ( 2.46 - 2.47 ), we get 
the following expression for the [C^^] matrix. 



_ [C]aa^[C] 

^ A + aT^[C]a 


{l-D). 


(3.25) 


3.3.2 Updated Lagrangian Formulation for Geometric Non- 
linearity 

When a body experiences large deformation and/or rotation, the equilibrium must be 
established in the current configuration. The normal solution procedure adopted in this 
case is updated Lagrangian. That is to say that when all the kinematical variables axe 
known from time 0 to time t in discreate time steps, the objective is to establish the 
equilibrium in the configuration at time t At . 

In order to derive the finite element equilibrium equations at time t + At, the 
principle of virtual work requires that 

^ (3.26) 

where is the Cauchy stress tensor at time t + At, is the infinitesimal 

strain tensor at time t + At, is the virtual work of the external body and surface 

forces at time t + At, is the domain at time t-\r At and is a volume element 
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at time t + At. The major difficulty in application of equation ( 3.26 ) is that the 
configuration at t + At is unknown. Moreover quantities like Cauchy stress tensor are 
not purely additive as now one has to take care of rotation into account. 


An elegant way of formulating the problem is given by Bathe [ 47 ]. Here, the virtual 
work expression ( equation 3.26 ) is written in terms of 2nd Piola Kirchoff stress tensor 
and Green Lagrange strain tensor which are energetically conjugate to each other. Thus, 
for updated Lagrangian formulation, the virtual work expression is written as 



S dv^ = 


(3.27) 


The right superscript stands for the current configuration and the left subscript stands for 
the reference configuration. Sij denotes 2nd Piola Kirchoff stress tensor whose relation 
with Cauchy stress tensor cT,j is given by 


Ct+At _ P _i+At t 


(3.28) 


Since Cauchy stress tensor is always referred to the current configuration, the left sub- 
script has been omitted. Here represents the derivative of the position vector 

at time t with respect to one at time t + At and p^l denotes the ratio of densities 
at time t and t + At. The Green-Lagrange strain tensor is defined as 


„i+At 


5( 




+ 


..t+Ai 


+ 


j.tH-At 

t^k,i 


i^k,j ) 


(3.29) 


where is the derivative of the displacement vector from time t to t + At with 

respect to the position vector at time t. 


An incremental decomposition of stress and strain gives 


tOij 




= orb -l-t A5,j, 

(3.30) 

J+At 

— t Acij "l-t At^ij, 

(3.31) 
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where, 


iA.c,j — 2 "ht ^Uj^i) 


(3.32) 


and, 


tArjij = - tAuk,i tAukj. 


(3..33) 


Here tAu,-j stands for the derivative of tAu (incremental displacement vector at time 
t) with respect to tx (position vector at time t). Thus, equation( 3.27 ) can be written 
with incremental decomposition as 


tASij 6{ tAcii) dv^ + tASii 6{ tArjij) dv^ + ^ 6{ tArjij) dv^ = 

(3.-34) 

The 2nd integral on the left hand side is a higher order term compared to other terms 
and hence can be neglected. More over, using an elastic-plastic incremental stress-strain 
relationship (equation 3.24), the 1st term can be approximated as .C.fS .Aeii 
The left subscript of denotes the time at which it is to be evaluated. The above sim- 
plification will definitely give rise to an error in the right hand side of the equation ( 3.34 ). 
The error is normally neutralised by using some iterative scheme like Newton Raphson 
Method. Thus, equation ( 3.34 ) can be written as 

tCfjki tAcki S{ tAtij ) du* -t- S[ tArjij ) du* = 

4 6{ tAcii ) dvK (3.35) 

Owing to the symmetries in C^kh 7«i equation ( 3.35 ) can be rewritten 

in the following form; 

(S t Ae^) tAe du* -f (S t A4) [TY t A^ du‘ = 

~ j (S j Ae^) a! du‘ 


(3.36) 
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where, 



tAe22^ tAe33, 2 

tAei2i2 ; 

t Ae235 2 

t A€3 i)^, 

(3.37) 

2:‘ = 

/ _t A 

V^ll? ^22 

5 <^33^ 

<^125 

^23? 

\T 

‘^3l) 5 


(3.38) 

tAur 

3 t A 1^24 

f Au2. 

,2 tAu2,3 

i Au3,i 

tAuz,2 tAu.3,3 

f. (3.39) 



■ 

0 

0 ■ 





[rr = 

0 

E‘ 

0 



(3.40) 


and 


E = 


L 0 

0 

r J 


<^12 


to 

<^22 

^23 

L ^31 

^32 



(3.41) 


The matrix is the elastic-plastic matrix evaluated at time t which is given by 

equation ( 3.21 ) or (3.25). 


3.4 FEM Equations 


The domain is discretised into a number of elements and the displacement field within 
each element is approximated as 


j Au = [Q] f Au® 


(3.42) 


where, 

tAu® = (AuJ, Auj) Aug, ... ,Au”, Am^j 


and Alt j, Auj, Aug stand for incremental displacements of node ’i’ in a;, y, s directions 
respectively. Here 




(3.43) 
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where, 


^ = [Nx 0 0 N 2 0 0 ... Nn 0 0], 

= [0 0 0 A^2 0 . . . 0 iV„ 0], 

^ = [0 0 ATi 0 0 7 V 2 . . . 0 0 


With displacement field defined by equation ( 3.42 ), the strain field within the 
element can be calculated in terms of nodal displacements. From equation ( 3.42 ), 


where 


d t^Ui IP ^ 

t^Uij — _ — tQ. t^u 

a tXj — 




(3.44) 


(3.4.5) 


'’J d tXj ' 

Substituting equation ( 3.44 ) in equations ( 3.32 ) and ( 3.33 ) and arranging them in 
vector form, one can get the strain displacement relations; 


where, 


and 


jAe = t[B]L tAyf, 
tAr? = t[B]N tAu^ 


t[B]L 




tQ 


T 

.2,2 


‘^2,1 


(3.46) 

(3.47) 


(3.48) 


2 


(3.49) 


Using the strain displax:ement relations derived above, equation ( 3.36 ) can be 
expressed as 
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AB]l AC?” AB]l dv^ ,Au4 
+ E{«(.A!i')’' .1^15 [T]‘ ,(£]„ dv‘^ ,Aii'| 

= J2[d{,Au’fr -SUd^u^f AB]1 z' (3.50) 

Here the term has been expressed in terms of the elemental nodal forces Ff using 

a standard procedure. is the elemental domain at time t. Expressing the term within 
1st parenthesis as linear elemental stiffness [ and that of 2nd one as nonlinear 

elemental stiffness [ Y above equation can be written as 

5 ; ([A'i]' + [AVd') ,Aji' = - /') (3.51) 

e e 

or 

i:[AT,A!,'=5;(£'-r) (3-52) 

e e 

where 

r = X, AB]lz‘dv‘ (3.53) 

gives the internal force vector. Assembling the elemental matrices [A']® and the elemental 
vectors F® and /® over all the elements, we get the following equation. 

[F],Au = (£-/). (.3.54) 

Here [K] is the global stiffness matrix, tAu is the global incremental displacement vector 
( at time t ) and F and / are the global external and internal force vectors respectively. 

Equation( 3.54 ), which represents the equilibrium equation in the current configu- 
ration, has to be solved iteratively. In an iteration cycle. The R.H.S of equation ( 3.54 ) 
represents the unbalanced force. The iterations are continued till the R.H.S becomes zero. 
The details of various iteration schemes are described in section 5. In each iteration, first 
the incremental strain tensor tAe,j is calculated from the strain displacement relation 
(equation 3.46 ). Then the Jaumann stress rate is calculated from equation (3.24) as 

t^€ki. (3.55) 
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Finally the Cauchy stress tensor at t + At is evaluated from 

( 3 - 56 ) 

where 0,]^ is a spin tensor defined by 

= tAui.,). (3.57) 

3.5 Numerical Scheme 

The finite element equation of an elastic-plastic material (equation 3.54 ) can be solved 
iteratively by various numerical schemes. Newton Raphson scheme is one of the most 
commonly used method for such problems. However, if the incremental load step is too 
large, the error involved in the computation of stress trom the matrix, evaluated 

at the previous point, also becomes large [48]. In such a situation, the backward return 
algorithm can be used to minimize the error. Another limitation of Newton Raphson 
method is that it often fails to converge in the neighborhood of critical point. In such a 
situation, the arc length method can be used. These three methods ( Newton Raphson 
method, the backward return algorithm and the arc length method ) have been discussed 
in the next three sub-sections. 

Normally we encounter two types of problems in mechanics: (1) Load control prob- 
lem and (2) Displacement control problem. In a load control problem, the desired defor- 
mation is achieved by prescribing a load at a point. On the other hand, in a displacement 
control problem, it is the prescribed displacement which gives the desired deformation. 
In a load control problem, if the load falls due to change in geometry, then most iterative 
schemes fail to converge. Displacement control problems are free from this limitation. In 
this thesis, only displacement control problems are considered. 


centrh umfkm 


I* 
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3.5.1 The Newton Raphson scheme 

The scheme is best described with the help of Fig { 3.1 ) in the following steps. In Fig 
( 3.1 ), m is an equilibrium point at load level P . The aim is to reach another equilibrium 
point (m + 1) at a load level (P + AP). 

Stepl— First of all, the stiffness matrix [A'] in equation ( 3.54 ) is evaluated using 
the stresses corresponding to point m. Next, equation ( 3.54 ) is solved with / = 0 and 
P = AP to get an incremental displacement vector t Au. 

Step 2- From the incremental nodal displacements, incremental strain components 
are calculated by using equation (3.46 ). Incremental stresses are calculated from the 
incremental strains by using suitable constitutive relation. Generally incremental stresses 
are calculated at the Gauss points of the element. If the point under consideration is in 
elastic range, the elastic stress strain matrix [C] is used. Otherwise, elastic-plastic matrix 
[C^^] is used. Since elastic-plastic matrix is dependent on the current state of stress, 
the stress value at position m is used for the computation of elastic-plastic matrix. The 
error involved due to this approximation is neutralised by backward return algorithm 
discussed later. Equation ( 3.55 ) is used to calculate Jaumann stress rate and the total 
Cauchy stress follows from equation ( 3.56 ). 

Step 3- Once the stresses are calculated at all the Gauss points of all the elements, 
the internal nodal force vector / is calculated. This force corresponds to the equilibrium 
load at point i in Fig( 3.1 ). 

Step 4- Since the applied nodal forces (P + AP) are not equal to internal nodal 
forces /, the deformed configuration is not an equilibrium configuration at load level 
(P + AP). It gives rise to an unbalanced nodal forces (P + AP — /). 

Step 5- The steps 1 to 3 axe repeated with unbalanced nodal forces as the load 
vector and i aa a starting point. Successive repetition of the above steps will lead to the 
equilibrium point (m-fl) for load level P-f AP when the unbalanced nodal forces become 
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almost zero (less than a pre assinged quantity). In each successive iteration, incremental 
displacement vector is added to its previous value. 

Faster convergence is achieved if the tangential stiffness matri.v at each successive 
equilibrium points is updated. The general practice is to achieve the total load increment 
in several steps. Within each incremental load step. Newton Raphson iterative scheme is 
applied with the same tangential stiffness matrix. Stiffness matrix is updated only when 
going for the next load increment. This modified Newton Raphson scheme is a good 
compromise between reduced computational time and faster convergence. 

3.5.2 The Backward Return Algorithm 

Consider the increment from point m to point i. In calculation of the stresses at point i 
( in step 2 of the Newton Raphson scheme ), we use the incremental stress-strain relation 
( equation 3.4 ) where the [C^^] matrix is evaluated at point m. When the incremental 
load step is too large, this gives stresses not at point i but at some point i which is not 
on the yield surface ( see Fig 3.2 ). Thus 

SLi' = £m + Ae. (3.58) 

A further return is required to bring the stresses to point i. Therefore we can write 

£.i — + Am (3.59) 

To determine Act, we proceed as follows. From equations (3.18) and ( 3.21 ) we note that 

= [^] - ^rndt [C] a,„. (3.60) 

Therefore, to bring the stresses to point i, a correction is needed only in the second term. 
We assume that the correction can be approximated as 

Ag:=-Xodt[C]a,>. (3.61) 

To deternoine the scalar Aq, we expand the plastic potential at point i . Thus. 
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Analogous to equation (3.12), we define A-- as 


Ay = : — ~Ap. 


(3.63) 


Xodt dp 

Substituting equations (3.61) and (3.63) in (3.62) and using the definition of flow vector 
( equation 3.11 ), we get 

Fi = F-> — Xodt gJi [C\ g^' — Xodt A-/. (3.64) 

Since the plastic potential has a zero value on the yield surface, F,- = 0. Thus 

F;- 


Xodt 


gj, [C\ g-i + .4;' 


(3.65) 


Substitution of equations (3.61) and (3.65) in (3.59) gives the following approximate 
expression for the stress at point i: 




Fj’ [C| a/ 

aj [C] Oj. + Af 


A further expansion will lead to the result 

F. [C] a,. 




F,» [C] a,. 


(3.66) 


(3.67) 


^ [C*] £Lj-' + oT/ [C\ g^ii + A,-" 

where the point i" is closer to point i than i' ( see Fig 3.2). When the plastic potential 
at point i becomes less than a pre assigned small number, no additional terms need to 
be added to equation ( 3.67). 


As a demonstration. Fig ( 3.3 ) shows the true stress versus true strain curve for 
a cylindrical specimen with and without backward return algorithm. The true stress 
is calculated by dividing the load with the current cross sectional area at the necked 
region. The true strain ( logrithmic strain ) is calculated from the formula 

e = 2\n(^j^ (3.68) 

where do and d are the initial and current diameters of the test specimen in the necked 
portion. Comparison of the curves with the actual curve shows improvement with the 
backward return algorithm. 
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3.5.3 The Arc Length Method 

The arc length method [ 49 ] is described here for a displacement controlled problem in 
conjunction with Newton Raphson method in its standard or modified version. Propor- 
tional loading is assumed but it can be extended to non proportional loading with minor 
changes. 

Notation : 

A left superscript indicates the configuration. Therefore, "^u, ....etc denote the total 

displacement vector, the load vector etc at configuration m. For proportional loading, 

the loads may be expressed by a single load factor ™A. Thus. 

™P = ^XP (3.69) 

where P is a basic load vector. Within one increment from configuration m. to (rn + 1), 
the positions i and j = (z -f 1) indicate the beginning and the end of an iteration cycle 
as shown in Fig ( 3.4 ). The changes in quantities from i to j ( both iterative points ) 
are denoted by and AA^l respectively. 


Starting point : 

The incremental equation ( equation 3.54 ) at position i may be expressed as 



‘[/i] = AP<^’^ -f ‘P- 7. 

(3.70) 

Denoting the out of balance force by i.e. 



'R= ‘£- 7 

(3.71) 

and using 

AP^^'> = AA<j7 

(3.72) 

we get 

^[K] Am^'^ = AA(^)P + ‘P- 

(3.73) 


Procedure for Displacement Control 

A simplified method as suggested by Batoz and Dhatt [ 50 ] is described here. Since the 
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equation ( 3.73 ) is linear, one can decompose the displacement vector into two 

parts 35 follows: 

(3.74) 

where 

‘[K] Au(^')^ = P (3.75) 

and 


‘[K] = ^R. 


(3.76) 


The nth component of the incremental displacement vector AuJ,^^ can be expressed as 

Au|f’’ = AA(J)Au|/')'^ + (3.77) 


Therefore, 


AA^^'^ 


A'uf/^ — 

n n 


(3.78) 


At the start of the iteration, every thing corresponds to the equilibrium point m. 
Therefore, °R — 0, = 0 ( from equation 3.76 ) and thus, = 0. Hence, 

AulH 


A(^) = AA(^) = 


(3.79) 


If we chose the nth component of the displacement vector as the prescribed displacement 
for the increment then, A^^^ becomes the factor for scaling the prescribed displacement 
from the displacement corresponding to the choice of basic load vector P_. 

For other iterations, = 0 , since is unaltered during iteration cycles. Therefore, 

AuW^ 


AA(^) = 


Aun 


U)I 


(3.80) 


Note that Au^-’^ and can be calculated by solving equations (3.75 ) and 

(3.76 ). Then, knowing AA^-'^ from equation ( 3.80), one can calculate the displacement 
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increment with the help of equation (3.74 ). Then the total clislacemeut vector at 

iteration point j becomes 

'u + Au^^K (3.81) 

Further, the value of A at point j is given by 

^A= b\ + AA(^>. (3.82) 

From this one can calculate the total load vector at point j as 

'P= 'AP. (3.83) 


The stress, strain and other secondary quantities can be calculated from this data using 
the standard relations. 
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Fig. 3.1 Illustration of Newton Raplison iterative scheme 
for a single load. 
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Fig. 3.2 : Illustration of Backward return algorithm 
for 1-D problem. 
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Chapter 4 


CASE STUDIES 


4.1 Introduction 

The updated Lagrangian finite element formulation of a damaged elastic-plastic material 
developed in the previous Chapter has been applied to a number of cases. The objective 
of these case studies is to investigate the effects of triaxiallity - plastic strain and 
damage variable on ductile fracture process. Finally, using the micro crack initiation cri- 
terion based on Thomason’s limit load model ( Chapter 2 ), a failure curve is constructed 
and a ductile fracture parameter is identified which is independent of geometiy. 


4.2 Finite Element Implementation 

The FE formulation developed in Chapter 3 has been implemented in an iterative fashion. 
That is, equation ( 3.54 ) is solved for equilibrium till the unbalanced force in anj’ iteration 
cycle is zero. The Arc length method in conjunction with Newton Raphson scheme 
( section 3.5 ) is used for displacement controlled problems. 

A step is defined as the collection of increments of all the quantities between two 
successive equilibrium states. W^ithin a step there are several iteration cycles. At each 
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step, the stiffness matrix is calculated by adding its linear and non linear components 
( equation 3.51 ). Within a step the stiffness matrix is kept constant and is updated 
only at the end of each step. Although the stiffness matrix is kept constant for iteration 
cycles, the internal nodal forces ( equation 3.53 ) and hence the unbalanced forces for each 
element axe calculated with a [B\i matrix evaluated at the current configuration. For this 
purpose, nodal coordinates are updated by adding incremental displacement components 
to achieve the current configuration. Following procedure is adopted for a iteration cy- 
cle. Incremental strain components are calculated from incremental displacements using 
equation ( 3.46 ). Next the Jaumann’s stress rate is calculated from equation ( 3.55 ). 
The total Cauchy stress follows from equation ( 3.56 ). The total displacement increment 
for the step is achieved by adding the incremental displacements for the iteration cycles. 

At the end of each step, the following quantities are calculated. The incremental 
strain tensor and the incremental equivalent strain for the total step are calculated from 
the total displacement increment for the step. The equivalent strain, thus obtained, is 
added to its accumulated value in the previous step to get the total equivalent strain. 
As the elastic strain is small, the total equivalent strain is treated as the total equivalent 
plastic strain. The equivalent stress (cTe,), the mean stress {(Jm) and hence the triaxiality 
(^) are calculated from Cauchy stress components. Incremental damage is calculated 
at the end of the step using the damage growth law ( equation 2.60 ). The incremental 
damage so obtained is added to its previous value to get the total damage. During 
elastic unloading, the calculations regarding equivalent plastic strain and the damage are 
bypassed. For reloading, the calculations for equivalent plastic strain and the damage are 
restarted from the point of unloading. The load is calculated from the nodal reactions. 
The programme is terminated when the condition for micro crack initiation by void 
coalescence (equation 2.80 ) is satisfied. 
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4.2.1 Validation of the Programme and the Convergence Study 

The elastic-plastic FE programme developed in Chapter 3 is validated with TEPSAC 
code developed by Tsu [ 51 ] for a test problem given in ref [ 51 ] . Fig ( 4.1 ) shows the 
geometry, boundary conditions and material properties of the problem. Fig ( 4.2 ) shows 
the load vs radial contraction at the mid-section of the specimen. The FE result obtained 
from our code matches very well with that given in ref [ 51 ]. The Maximum discrepancy 
is only 3 % which occurs at the maximum load. It is not necessary to perform the damage 
calculations during i.his validation study. It should be noted that the stress-strain curve 
given in ref [ 51 ] is approximated by the relation ay — 

The convergence study is done for the problem of a cracked plate in plane strain 
condition. Fig ( 4.3 ) shows the geometry and the loading of the problem. The load 
point displacement and the crack opening displacement (COD) at a point 0.2 mm away 
from the crack tip, are checked for different crack tip element sizes at the same load 
level. Table 4.1 shows the results of the convergence study. Fig ( 4.4 ) show its graphical 
representation. It is found from the convergence study that for a crack tip element size 
not exceeding 0.2 mm, the results are satisfactory in engineering sense. 


4.3 Case Studies 

4.3.1 General Description of the Problems 

The material used in this study is AISI-1090 spherodised steel whose void growth curve 
for a cylindrical specimen in uniaxial loading is given by LeRoy et al [ 36 ]. For the sake 
of completeness it is reproduced here in Fig (4.5 ). The chemical composition and the 
mechanical properties of the material are given in Tables 4.2 and 4.3 respectively. The 
coefficients oi , 03 and c of the damage growth law ( equation 2.60 ) are obtained as 
follows. First, equation ( 2.54 ) is rearranged as a relation between and e% using 
Bridgman’s [ 52 ] expression for triaxial:ity as a function of strain. Then, from 
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experimental results of LeRoy et al (Fig 4.5 ), the slopes are calculated at different 
strain levels. Finally, the coefficients ai , a 2 and c are obtained by the method of least 
square curve fitting. The following are the values obtained. 
ai = 9.8 X 10-°^ MPa-\ 

02 = 1.84 MPa-\ 

c = 1.898 X 10-“. (4.1) 


Five specimens with different geometries are studied. The geometries of the speci- 
mens are shown in Fig ( 4.6 ). For all cases except for the cracked plate, a total of 105 
eight noded isoparametric elements are used. For the cracked plate, the number of eight 
noded isoparametric elements is 182 including the crack tip elements. The size of crack 
tip elements is 0.2 mm x 0.2 mm. The deformed and the undeformed mesh patterns for 
the cylindrical specimen are shown in Fig ( 4.7 ). It is worth noting here that a small 
imperfection of the order of 0.1% is introduced in the diameter of the mid-section to 
simulate necking [53] (see the Appendix). 

4.3.2 Results and Discussion 

To check the accuracy of the result obtained from computer simulation, the simulated 
results for the cylindrical specimen are compared with the e.Kperimental results of LeRoy 
et al [ 36 ]. Fig ( 4.8 ) shows the simulated as well as experimental damage growth curve. 
The simulated curve is obtained by integrating the damage growth rate law up to the 
failure point. Note that the simulated curve depicts the experimental trend quite well. 
Further, it is able to predict the damage values with a reasonable engineering accuracy. 
Fig ( 4.9 ) shows the comparison between the simulated true stress-true strain curve 
with the experimental results of LeRoy et al [ 36 ] for the cylindrical specimen. The 
true stress is calculated by dividing the load with the current area of the necked region. 
The true strain is calculated from equation (3.68). It is observed from Fig ( 4.9) that 
the simulated stress-strain curve matches with the experimental curve very well ( The 
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difference in stress level is only 8% even at a strain level as high as 55% ). This shows 
that the proposed damage growth law (equation 2.60 ) is consistent with the experimental 
material behavior. 

Fig ( 4.10 ) shows the damage growth with equivalent plastic strain up to failure 
( i.e. up to micro crack initiation ) for the five cases. .A.s expected, the damage grows with 
plastic strain. For all cases except for the cracked plate, the maximum damage occurs at 
the centre. For the cracked plate, it occurs at the crack tip. The critical values of damage 
(Dc) and strain (e^^)c are obtained when equation ( 2.80 ) is satisfied. The values of the 
critical strain and the critical damage for the cylindrical specimen are respectively 60.6% 
and 5.62%. The reported value of the fracture strain for the cylindrical specimen in ref 
[ 36 ] is 63%. Since the critical strain is close to the experimental value of fracture strain, 
it shows that the Thomason’s model is a good representation of micro crack initiation 
due to void coalescence. 

Fig ( 4.11 ) shows the damage versus triaxiaLity curves for the five cases. It is 
found that, for the major portion of the deformation, triaxial jty remains constant. It 
means the proportional loading condition has been achieved in all the five cases. This 
situation arises because the loading is uniaxial tension in all the five cases. Since, the 
triaxiality function /(f^) does not change in damage growth rate equation ( 2.60 ), the 
dialatational rate becomes constant. Hence, the nonlinearity in damage growth relation 
is due to plastic strain only. Although the triaxial Ity remains constant in each case, its 
value differs from case to case since it saturates to different critical values for different 
cases. The maocimum critical triaxial ity is found at the crack tip of the cracked plate 
whereas the minimum value occurs for the cylindrical specimen. For the cracked plate, 
triaxial?ity increases in the initial part of the deformation due to change in crack tip 
geometry. 

Table 4.4 lists the critical values of strain, triaxial ity and damage for the five dif- 
ferent cases. It is found that for a wide variation of critical strain and triaxiallity, the 
value of critical damage remains almost constant. Hence, the critical value of damage 
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{Dc) can be regarded as a geometry independent parameter which indicates the micro 
crax:k initiation in ductile solids. 

Fig(4.12 ) shows the variation of triaxiality with equivalent plastic strain for all 
the five cases. The locus of the failure (or critical) points is called as the failure curve. 
The failure curve has a qualitative similarity with that of Hancock and Brown [ 54 ] but 
the quantitative comparison is not possible since they have used a different steel. This 
curve shows that the failure can not be predicted either by critical strain alone or by 
critical triaodalj ty alone. Since one has to take both into account, it is better to predict 
it by critical damage which incorporates both strain and triaxial tty. It is observed that 
critical strain increases with decreasing critical triaxiah ty. Finally, Fig ( 4.12 ) shows 
that, besides strain, triaxiall’ty also plays an important role in void growth and micro 

crack initiation. % 

1 


■The content of Chapter 2, 3, and 4 has been accepted in Engn. Fract. 
Mech. for publication. 
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Material : Properties 


Youngs modulus, E =71724 MPa 
Poisson ratio, ^ = 0.33 
Initial yield strength, 0^ =350 MPa 
Hardening parameter, K= 52370 MPa 
Hardening exponent , n = 0.07229 


Fig. 4.1 Geometry and boundary conditions of the 
test problem . 



Radial contraction (mm) 


030 


Fig 4.2. load versus radial contraction for tlie 
test problem 
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P 





All dimensions are in mm. 

Material ! AISI 1090 spherodised steel. 
Youngs modulus E = 210 GPa 
Poisson ratio,V=0.3 
Yield stress 0^= 464 MPa 
Hardening parameter K=1115 MPa 
Hardening exponent, n = 0.2. 


Fig. 4.3 Geometry and loading of the (1/2 TCT specimen 
for convergence study. 
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Fig 4.4. Graphical representation of conver- 
gence results. 
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Fig 4.5. Experimental void growth curve, ref [ 36 ]. 
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c 

Mn 

P 

S 

Si 

Fe 




0.022 

0.20 

balance 


Table 4.2: Chemical composition (% weight) of Material AISI-1090 steel 


poisson’s 

Young 

Yield 

Ultimate 

Fracture 



ratio 

modulus 

stress 

stress 

strain 



P 

E 

<^0 

^uts 


K 

n 


GPa 

MPa 

MPa 


MPa 


0.3 

210 

464 

t 

619 

0.63 

1115 

0.19 


Table 4.3: Mechanical properties of Material AISI-1090 spherodised steel 

















(a) Cylindrical 

(b) Cylindrical prenecked 

(c) Plane strain 

(d) Plane strain side groovec 

(e) Cracked plate 

All dimensions are in mm 



Fig. 4.6 specimen geometry 






28 mm 



Fig 4.7. Deformed and imdeformed mesh pat- 
terns for cylindrical specimen. 
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Plastic strain 


Fig 4.8. Comparison of computer simulated 
damage growth curve with experimental re- 
sults for the cyhndrical specimen. 




Fig 4.9. Comparison of computer simulated 
true stress-true strain curve with experimen- 
tal results for the cyHndrical specimen. 
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True stress ( MPa ) 



1500 







Damage (^) 


I Cylindrical 2 Plane strain 

3 Cylindrical prenecked 4. Plane strain side grooved 
5. Cracked plate 


Triaxial ity (^) 

pn 


Damage versus triaxial ity curves for the 
five specimens. 
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SI No 

Specimen 

Critical strain 

Critical triaxiaUty 

Critical damage 

1 

cylindrical 

0.61 

0.593 

5.62 X 10-°2 

2 

plane strain 

0.55 

0.701 

5.68 X 10-°2 

3 

cylindrical prenecked 

0.50 

0.820 

5.61 X 10 -°^ 

4 

plane strain side grooved 

0.46 

0.925 

5.90 X 10-°2 

5 

cracked plate 

0.34 

1.430 

5.80 X 10-°2 


Table 4.4: Critical values of quantities for the five different cases 
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1. Cylindrical 

2. Plane strain 

3. Cylindrical prenecked 
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Plastic strain 


Fig 4.12. Failure curve. 



Chapter 5 


APPLICATION TO FRACTURE 
MECHANICS 


5.1 Introduction 

The previous Chapter was devoted to identify a continuum parameter namely damage 
variable {D) to indicate micro crack initiation in ductile material. It has been shown that 
the value of the damage variable reaches a critical value at the micro crack initiation. 
It has also been pointed out that, for different geometries, even though there is a wide 
variation in the values of critical strain and triaxial ity, the value of the critical damage 
remains almost the same. 

In this chapter, a local criterion for crack growth initiation ( the onset of growth of a 
existing crack ) in Mode-I ductile fracture is proposed using the critical value of the dam- 
age variable (Dc). To test the validity of this criterion, experiments have been conducted 
on standard specimens and have also been simulated numerically. The proposed criterion 
has been used to find the critical values of load (Per), fracture toughness ( J/c) and crack 
tip opening displacement ( iCTOD)c )• The critical values of the fracture parameters, 
obtained from computer simulation, compare favourably with their experimental values. 
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5.2 Crack Growth Initiation Criterion 

A local fracture criteria can be stated by two parameters namely a critical continuum 
parameter and a critical length parameter. The critical length parameter denotes the 
distance ahead of the crack tip where the continuum parameter has to reach its critical 
value for micro crack initiation. Thus, the critical length is a characteristic dimension 
of the material which fractures when the crack growth initiates. It is obvious that this 
characteristic dimension depends on the microstructure of the material and hence it is a 
material property. 

In Chapter 1, some local crack growth initiation criteria have been discussed. In the 
RKR ( Ritchie, Knott and Rice [42] ) criterion, stress is used as the continuum parameter 
while the critical length parameter is taken to be twice the average diameter of Fe grain. 
On the other hand, in the criterion of Ritchie, Server and Wullaert [43 ], strain is used 
as the continuum parameter. The critical length parameter is coming out to be one to 
six times the average spacing of MnS particles. It is possible that in a few cases, the 
critical stress or critical strain may be able to predict micro crack initiation in ductile 
materials. However, as stated in Chapter 4, in general, micro crack initiation depends not 
just on the critical strain or critical triaxiality but on a critical combination ot these two 
parameters. Since the damage incorporates both these parameters, in the present study, 
the critical damage {Dc) is used as the critical continuum parameter. For the material 
AISI-1090 spherodised steel, the value of Dc for different specimens lies between 0.056 
to 0.059 ( see Table 4.4 ). In further studies, this value is taken as 0.05. In the present 
study, the critical length parameter (Ic) is taken as the average austenite grain size of 
the material. The choice of austenite grain size is not unjustified because the larger sized 
particles restrict the austenite grain growth. Thus the proposed crack growth initiation 
criterion can be stated as follows. Whenever 

D = Dc a.tl = lc (^-1) 


the crack growth initiates. 
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5.3 Experiments 

The experiments were performed on AISI-1095 spherodised steel whose chemical compo- 
sition is given in Table 5.1. The specimens tested for experimentation are Three Point 
Bend ( TPB ) specimens and Compact Tension ( CT ) specimens. The dimensions of the 
specimens were chosen according to the ASTM standard as far as possible. Figs ( 5.1(a) ) 
and ( 5.1(b) ) show the geometry and dimensions of the specimens. 

The bulk material was austenised by heating it to 1000°C for two hours and then 
furnace cooled. Fig ( 5.2 ) shows the austenised grain structure of the material. The 
average grain size was measured by interception method. The average grain size was 
found as 33/xm within a range of SOfim to Next, the material was supplied for 

the preparation of specimens. After machining, the specimens were annealed by heating 
them to 750°C for two hours to relieve the residual stresses due to machining. The final 
heat treatment was done to spherodise the material by heating at 710°C for 20 hours. 
Fig ( 5.3 ) shows the spherodised structure of the material. 

A tension test was performed to find the mechanical properties of the material. 
Table 5.2 shows the results of the test. 

The specimens were precracked by fatigue loading in 100 KN MTS machine. The 
maximum load for a fatigue cycle was kept within 60% of the limit load which was 
calculated from the following formulae [ 55 ]. 

for TPB specimen (5-2) 

3 S 

and, 

p — for CT specimen. (5-3) 

^ 2(W-ha) 

Here, S is the half-length of the TPB specimen, B is the width of the specimen, b is the 
ligament length, a is the effective crack length and = a -f 6 { see Fig 5.1 ). After the 
appearance of a crack, the maximum load was reduced to 40% of the limit load gradually. 
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Following are the fatigue data for a TPB and a CT specimen. 

For a TPB specimen, 

Initial ma.ximum load = 8 KN, 

Fatigue ratio = 0.1, 

Frequency = 10 Hz, 

Final maximum load = 5 KN. 

A total of 1,10,000 cycles were required to have a fatigue crack growth of 3.5 mm (which 

was measured after breaking the specimen at the end of the experiment). 

For a CT specimen. 

Initial maximum load = 6 KN, 

Fatigue ratio = 0.1, 

Frequency = 10 Hz, 

Final maximum load = 4 KN. 

A total of 66,000 cycles were required to have a fatigue crack of size 2.3 mm . 


After precracking, the specimens were steadily loaded by controlling the displace- 
ment of the ram. In all the cases, a ram speed of 0.1 mm/min was maintained. The load 
ajid the load point displacement were recorded in X-Y plotter throughout the loading 
history. The crack growth initiation was indicated by a distinct pop-in in all the cases. 
The measured values of the load and the load point displacement and the calculated value 
of J/ at this point are termed as their critical values and are denoted respectively by P^., 
Ucr and Jic • For calculating J/ value at a certain load level, the following formulae were 


used [56 ]. 

For TPB specimen. 



(5.4) 


where A is the axea under the load versus load point deflection curve. For CT specimen. 


Ji = 


t}A 

Sb 


(5.5) 


where Ernest f<ictor r] is given by [56] 


r} = 2.0 + 0.522 


W 
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Figs ( 5.4(a), (b), (c) ) and ( 5.5(a), (b) ) show the experimental load deflection 
curves for the TPB and the CT specimen respectively. The fractured surfaces of broken 
specimens were scanned by 15 KV JEOL JSM 840 A scanning electron microscope. Figs 
( 5.6(b), (c), (d), (e) ) show the fractography taken at different places. The locations of 
scanning were chosen in accordance with Fig 5.6(a). 


5.4 Numerical Simulation Of Fracture Tests 

The updated Lagrangian finite element formulation of a damaged elastic-plastic material, 
developed in Chapter 3, is used to numerically simulate the fracture mechanics tests 
conducted in this study. The damage growth law, used in computing the damage variable 
’Z) ’ at any deformation level, is the same as given by equation ( 2.60 ). Since the chemical 
composition and the mechanical properties of the test material are close to that of AISI- 
1090 spherodised steel ( see Tables 4.2, 4.3, 5.1, and 5.2 ), the values of the coefficients 
ai, 02 and c of the damage growth law are taken to be the same as that of AISI-1090 
spherodised steel. These values axe given by equation (4.1). Also, the critical value 
of the damage variable is taken same as that of AISI -1090 spherodised steel i.e Dc = 
0.05. The average austenite grain size of the test material, as mentioned in the previous 
section, is 33^m. For crack growth initiation, the criterion proposed in section 5.2 is 
used. According to this criterion, whenever the value of the damage variable reaches 
Dc at a distance Ic ( austenite grain size ) ahead of the crack tip, the crack growth 
initiates ( equation 5.1 ). Size of the crack tip element is chosen so as to facilitate the 
application of this criterion. Fig ( 5.7 ) shows a representative crack tip element of size 
2/c X 2/c such that the Gauss point nearest to the crack tip represents one grain. Figs 
( 5.8(a), (b) ) show the finite element mesh for both TPB and CT specimens. Eight noded 
isoparametric elements are used everywhere except at the crack tip where the collapsed 
crack tip elements axe used. A 2 x 2 Gauss integration scheme is used throughout. 


The specimens are deformed by giving a prescribed incremental displacement at a 
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fixed point which is the load point in actual experiment. Load is calculated from nodal 
reactions. At the end of each increment, the increment in damage is calculated by using 
the damage growth law (equation 2.60) which is added to its previous value to get the 
total damage. During elastic unloading the damage calculation is bypassed. At the end 
of each increment, the damage at the Gauss point nearest to the crack tip is compared 
with its critical value. As soon as ’D’ reaches the programme is terminated. All 
values at this stage are termed as the critical values of the respective quantities. The 
’Jj’ value is calculated after computing the area under the load deflection curve using 
the relations (5.4 - 5.5). The stretch zone width is calculated from the stretching of the 
initial crack i.e 

SZW = (a,- - flo) (5.6) 

where, a,- is the current crack length and is the initial crack length. The value of 
the stretch zone width at the critical condition is regarded as critical stretch zone width 
or {SZW)c. Out of all collapsed nodes, the 1st node is vertically restrained and others 
are allowed to move during deformation (Fig(5.9)). CTOD is regarded as the vertical 
opening of the last collapsed node (node 9 ). Its value at the critical condition is regarded 
as the critical crack tip opening displacement or {CTOD)c- 

Out of all test specimens, TPBl and CTl were completely simulated using the 
procedure cited above. The plane strain condition was assumed in both the cases. 


5.5 Results and Discussion 

Experimental Results 

Experimental load deflection curves are shown in Figs ( 5.4(a), (b) 5 (c) ) for TPB 
specimens. It is observed that all the curves show a distinct pop-in which is taken as the 
crack growth initiation point. This is followed by a ductile crack growth. At the crack 
growth initiation, there is profuse void growth and coalescence in the highly localised 
plastic zone at the cradk tip and this causes the initial instability. This tendency reduces 
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and finally arrested due to fall in triaxiality as the free surface is approached. 

Figs ( 5. 5(a), (b)) show the load deflection curves for CT specimens. The trends are 
similar to that of TPB specimens. 

Table 5.3 shows the experimental results of TPB specimens. It is found that, for 
alw ranging from 0.526 to 0.595, J/c value varies from 65.77 KN/m to 70.99A'A^/m. The 
average experimental value of is 66.66KN jm. The Table shows a decrease in critical 
load with the increase of ajw ratio. This is due to decrease in the ligament length h. 
Table 5.4 shows the experimental results of CT specimens. The average value of J/c is 
found as 65.6SKN/m. 

The fractographs (Fig 5.6(b), (c)) reveal the critical stretch zone. The critical 
stretch zone width for TPBl specimen is found to be Z5fim and that of CTl specimen 
to be bOfxm. The other fractographs (Figs 5.6(d), (e)) show voids in the voided zone and 
ductile shear lips near the free surface. 

Comparison between Experimental and Numerical Results 

Fig (5.10) shows a comparison of the computer simulated load versus load point de- 
flection curve with that of the experimental curve for the TPBl specimen. The computer 
simulated curve depicts similar trend as that of the experimental curve. The difference 
in critical load prediction is 11.0%. Fig (5.11) is a comparison of load versus load point 
deflection curve with that of the experimental one for CTl specimen. The difference 
in critical load values is 5.0%. Both the computer simulated curves have higher critical 
loads than those of experimental curves because the simulated curves are for the plane 
strain condition. 

The computer simulated J/c value of TPBl specimen is 72.4AKNfm. On the other 
hand, the average experimental value of Jjc for TPB specimens is 68.08 KNjm ( Table 
5.3 ). The simulated value is 6.3% higher than the experimental value. The computer 
simulated value of J/c for CTl specimen is 73.78KNfm which is 11.0% higher than that 
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of the average experimental value (65.68A'iV/m) for CT specimens ( Table 5.4 ). Both 
the simulated cases yield higher values of than those from experiments. Since the 
computer simulations are for plane strain condition, the simulated values are expected 
to be less than those of experiments. The discrepancy may be due to the reason that the 
critical value of damage variable ( ) of the material is determined by using Thomason’s 

limit load analysis (Chapter 2 ) which uses the upper bound theorem. Thus, the value 
of Dc equal to 0.05 is an over estimation of the actual value. 

The values of critical stretch zone width, calculated for TPBl and CTl specimens 
from equation (5.6), are 31.26/im and 32.09^m respectively. These computer simulated 
values are almost equal to the average austenite grain size ( 33/tm ) of the material. The 
experimental values of critical stretch zone width ( 35/xm for TPBl specimen and 50/im 
for CTl specimen ) are within the range of measured grain size ( 30fj,m - 60/ur7z ). The 
results show that the critical stretch zone width is of the order of grain size Ic- 

Comparisons of computer simulated results with experimental data show that the 
overall agreement is satisfactory. Since the computer simulation uses the proposed crack 
growth initiation criterion ( equation 5.1 ) for calculating the critical values of the quanti- 
ties, it shows that the proposed criterion is reasonably good in explaining Mode I fracture 
for the class of materials discussed here. 

Additional Numerical Results 

Some additional quantities are computed from the plane strain analysis of a Full 
CT specimen (Fig ( 5.12 ) ) . Fig ( 5.13 ) shows the growth of plastic zone at different 
deformation levels. The critical plastic zone size is found to be 0.176 mm. The critical 
damage zone is shown in Fig ( 5.14 ). The ratio of the critical plastic zone to critical 
damage zone is found as 5.6. This shows that the damage is a highly localised event at 
the crack tip. Figs ( 5.15 and 5.16 ) show the variation of triaxiairty and (Tyy along the 
crack line. It is observed that both the values drop at the crack tip. This is expected 
due to crack tip blunting. 


9 


A graph of Ji versus CTOD is shown in Fig ( 5.17 ). The relation is linear upto crac 
growth initiation point. The value of critical crack tip opening displacement ( CTOD, 
IS 56/^m which is 1.8 times the calculated value of critical stretch zone width. Slope c 
the curve is found as ‘i.da,. The theoretical value of the slope is 1.6a, ( in plane strai. 
condition ) for a ideally plastic material in small scale yielding. In large scale yielding 
this value can go upto d.Ga, depending on hardening exponent (n) [57 ]. Thus, the valu 
of the slope is in confirmity with the exiting result. 
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c 

Mn 



S 



■ 

Si 

Fe 



i 




1 

1 


Table 5.1: Chemical composition (% wt ) of the test material 


Young’s 

Poisson’s 

Yield 

Ultimate 

Ultimate 

Fracture 

Hardness 

K 

n 

Modulus 

Ratio 

Stress 

Stress 

Strain 

Strain 




E 

1/ 

CTo 

O'uts 






GPa 


MPa 

MPa 



Rc 

MPa 


VJi Jl 

210 

0.3 

497 

690 

f 

15% 

59% 

34 

1211 

0.22 


Table 5.2: Mechanical properties of the test material 
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Three point bend (TPB) specimen All dimensions ai-e in mm 


Material 

Specimen 

No 

a 

mm 

b 

mm 

W 

mm 

B 

mm 

S 

mm 

AISI 1095 

Spherocliseci 

Steel 

TPBl 

11.9 

8.1 

20 

20 

40 

TPB2 

lO’O 

9.0 

19 

19 

40 

TPB3 

10.1 

8.9 

' 19 

19 

40 


Fig. 5.1 (a) : Geometry of TPB Specimens. 





Compact Tension (CT) Specimen All dimensions are in mni 


Material 

Specimen 

No 

a 

mm 

b 

mm 

W 

mm 

B 

mm 

AISI 1095 
Sphcroclisecl 
Steel 

CTl 

14.5 

10.5 

25 

8 

CT2 

15.0 

10.0 

25 



Fig 5.1 (b) : Geometry of CT Specimen^i. 





Fig 5.2. Aiisteiiite gr? 
AISI-1095 ( C% 0.92) 


Fig 5.3. Spherodised structure of material AISI- 
1095 steel ( C% 0.92). 
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Fig. 5.4 (a) Experimental load deflection curve of specimer 
TPBl. 
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0 0.5 1.0 1.5 2.0 


Load point deflection (mm) 

Fig. 5.4(b) Experimental load deflection curve of spec 
TPB-2.. 



Load point deflection (mm) 


Fig. 5.4(c) Experimental load deflection curve of sp 
TPB-3. 




0 0.5 1.0 1.5 2.0 2.5 


Load point deflection (mm) 

Fig. 5.5 (a) Experimental load deflection curve of 
specimen CT- 1 



Load point deflection (mm) 

Fig. 5.5(b) Experimental load deflection curve of 



Specimen 

Su 

w 

a, tsi 

Ucr 

mm 

Jlc 

KN/m 

Average J/^ 

KN/m 

TPBl 

0.595 

12.0 

0.65 

67.50 


TPB2 

0.526 

15.1 

0.60 

70.99 

68.08 

TPB3 

0.531 

15.0 

0.59 

65.77 


Table 5.3: Experimental results of TPI 

specimens 

Specimen 

w 


C/cr 

mm 

Jjc 

KN/m 

average 

KN/m 

CTl 

0.58 

6.98 

0.45 

64.52 

65.68 

CT2 

0.60 

7.70 

0.425 

66.85 



Table 5.4: Experimental results of CT specimens 






sz 

(a) Location of different zones. 



(b) Stretch zone of TPBl specimen. 



(c) Stretch zone of CTl specimen. 

Fig. 5.6 Fractographs of fractured surface. 





(d) Voided zone. 


(e) Ductile shear lips. 


Fig. 5-6 Practographs of fractured surface. 









Fig 5.8(b). Finite element mesh of specimen 

CTl. 





Fig. 5.9 Modelling of crock tip. 
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Fig 5.10. Load versus load point deflection 
curve for TPBl specimen. 
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Fig 5-11* Load versus load point deflection 
curve for CTl specimen. 





8.4 mm 








Chapter 6 


THREE DIMENSIONAL 
ANALYSIS 


6.1 Introduction 

In chapter 5, numerical estimation of fracture parameters has been done using the plane 
strain analysis. However, it is observed that the specimen thickness also plays an im- 
portant role in fracture process [ 58, 59, 60 ]. Truly speaking, fracture process is a three 
dimensional phenomenon as the deformation and stresses vary along the crack front. 
Even for a thick specimen, the plane strain condition prevails only in the middle part 
of the specimen whereas the surface regions remain in plane stress condition. Further, 
the crack front does not remain straight during propagation as it propagates more at the 
middle and less at the surface. Thus, both the crack initiation and propagation are highly 
dependent on the specimen thickness. In order to study the effect of specimen thickness 
on fracture parameters, a three dimensional analysis of the crack front deformation is 

necessary. 

The laboratory test data of fracture parameters pertain to only the plane strain con- 
dition. However, the plane strain and plane stress conditions are only approximations 
to the actual state of deformation/stress when certain conditions are broadly satisfied. 
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Thus, in vu'd.T !o apply fho laboratory test data to actual fracture 
a ki;.av!,-.!ur uf ihr.M‘ dimensional deformation/stress field ahead of 


mechanics problems, 
a crack front is nec- 


t'ssary. 


In th< « «ir!i« r Chapters, two dimensional ( mainly plane strain ) analysis of the 
ductile fracture process was carried out using the updated Lagrangian elastic-plastic 
forniuiatiun for a dantaged material. Here the same method has been extended to three 
disnensiiiiial anaiysis of the crack front. Due to limitation of computer resources, only a 
( rude me.sh is used for the three dimensional analysis. Therefore, the results obtained 
are of <|ualitative nature only. The objective of the study is to investigate the variation 
of the contimmin parameters along the thickness as well as the dependence of the critical 
fracture parameters on the thickness. 


6.2 The Finite Element Analysis 

1/2T CT specimens { Fig 6.1 ) with varying thickness ( B ) are analysed in FE study. The 
material used is A 1ST 1095 spherodised steel with chemical composition and mechanical 
properties given in Tables 5.1 and 5.2 respectively. Numerical schemes used are the same 
as that for the two dimensional analysis. Each specimen is divided in 3 layers across 
the half thickness. Each layer contains 84 twenty noded brick elements. The size of 
the elements along the crack front, in X-Y plane, is chosen to be of the order of the 
grain size of the material { see Section 5.4 ). Due to symmetry of the geometry and also 
of the loading condition, only a quarter of the specimen is analysed. Proper boundary 
conditions axe applied at the planes of symmetry. In the Z plane, the Z displacements 
axe restrained whereas in the Y plane, the Y displacements are restrained except at the 
crack face. In order to restrict the rigid body movement in X direction, a point A on the 
midplane which is far away from the crack front ( Fig 6.1 ) is restrained from moving in 
the X direction. A 2 x 2 x 2 Gauss integration scheme is used. All quantities ( stress, 
strain, triaxial ity and damage ) are calculated at the Gauss points. These values are 
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extrapolated to the mi.lplane while applying the crack growth initiation criterion. The 

. . ^ thjf^kness is calculated by summing the nodal reactions across the 

t in‘.s.s. I he load point displacement is taken as the average displacement over the 
11 riims, the value of .// calculated from the load deflection curve with the help 

ti. uiiuilivn I ) ,s an average value over the thickness. While calculating the crack 
tip mng duspLuenient, the opening of the crack faces at a distance of OMCwim ( one 
element size J behind the crack tip is considered. 


6.2.1 The Crack Growth Initiation Criterion 

crack growth initiation criterion proposed in section 5.2 for two dimensional case 
{ equation o.i ) }$ extended here to the three dimensional analysis. The midplane, 
perpendicular to the crack front, is chosen for application of the failure criterion. In 
this plane, D has to reach the critical value at the characteristics length /„ ( average 
austenite grain size ) ahead of the crack tip for the crack to grow further. Thus, for crack 
growth initiation 

D = Dc at / = 4 in the midplaiie. (6.1) 

Therefore, the dimensions of the elements along the crack front, in X-Y plane, are selected 
as 2la X 24 . The dimension of the elements in Z direction is chosen more at the centre 
and less at the surface to take care of the sharp variation of the stresses close to the 
surface. The values of fracture toughness ( Jj ), triaxiahity and strain at the moment 
when the above condition is satisfied are regarded as their critical values. 


6.3 Numerical Results and Discussion 

Fig (6.2) shows the comparison of load deflection curve obtained from 3-D analysis with 
experimentai and plane strain curves of Chapter 5 ( Fig 5.11 ). For the ease of comparison, 
these curves have been reproduced here. The specimen used in this study is a 1/2T CT 
specimen with thickness of 8mm. There is a close matching among the results. It means 
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thtit A pltijit* strain condition has almost been reached for a thickness of Smrn. The 
thickness requirement for plane strain condition, as recommended by ASTM. is 

S>25^, (6.2) 

where, i.s the plane strain fracture toughness and (Tq is the initial yield stress of the 
material. I his gives a required value of the thickness of the specimen as iSmm. Thus a 
thickness tjf bmin is well above the required value. 

Variation of Continuum Parameters along the Thickness 

V'arialion of certain continuum parameters along the thickness at the critical condi- 
tion of crack growth initiation is discussed here. 

Fig ( 6.3 ) shows the variation of CTOD along the crack front. It shows larger 
crack-tip opening at the mid-plane and less at the surface. 

Fig ( 6.4 ) gives the variation of triaxial-ity at the crack tip ( at a distance r ) 
along the thickness of the specimen. For the specimen of 8mm thickness, there is a 
small increase in triaxiality from the mid-plane towards the surface. The decrease in 
triaxiallity at the midplane is due to larger opening of the crack-tip at the mid- plane 
[59]. The triaxiality drops sharply near the free surface. For the specimen of 2mm 
thickness, triaxial ity falls gradually from the mid-plane to the free surface showing a 
smooth transition from plane strain condition to the plane stress condition. 

Fig ( 6.5 ) shows the variation of out of plane stress [cr^z) at the crack tip ( at a 
distance r ) along the thickness. For the specimen of 2mm thickness, is zero for 
a portion of the thickness near the surface indicating a complete plane stress condition 
near the surface. The other specimen of 8mm thickness has a high value of at the 
midplane with a fairly uniform distribution over a large portion of the thickness. The 
value of (Tzz falls sharply near the surface. From Figs ( 6.4 ) and ( 6.5 ), it appears that 
the state of stress in the specimen of 8mm thickness can be approximated by a plane 
strain condition whereas the specimen of 2mm thickness requires a three dimensional 
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analysis. 

t ig ( fi.f) I is a graphical representation of the variation of equivalent plastic strain 
at the rriw k tip ( at a distance r ) along the thickness. The thiner specimen (thickness 
of 2mm \ recpiires more strain for crack growth initiation than that of thicker specimen 
(tliickne.ss of -Smnj). The nature of the curves matches qualitatively with that given in 
ref I 58 1. 

The above discussion shows that both the distribution as well as the values of the 
continuum parameters depend on the thickness of the specimen. Thus, a ductile fracture 
process is a thickness dependent phenomenon. 

Criiiml Fraciurr Parameters 

Fig ( 6.7 ) shows the variation of critical Jr for specimens with different thickness. 
It shows that critical J; for 2mm thick specimen is almost 25% higher than that of the 
specimen of Smrn thickness. For a thidcness of 6mm and above, the critical J ; becomes 
almost constant and close to its plane strain value J/c- It should be noted that the plane 
strain condition is an ideal case. A fully plane strain condition can not be achieved even 
for a very thick specimen. However, an approximate plane strain analysis can be carried 
out for those specimen whose thickness satisfies the condition of equation ( 6.2 ). 

Fig (6.8) gives the variation of critical triaxial.tty ( in the midplane ) with the thick- 
ness. Thicker specimen shows larger critical triaxial ity. For the specimen of thickness 
6mm and above, the critical triaxial Ity saturates to its plane strain value. 

Variation of critical strain ( in the midplane ) with the thickness is shown m Fig ( 6.9 
). A thiner specimen requires higher critical strain. A 2mm thick specimen requires 16% 
higher strain than that of 8mm thick specimen for crack growth initiation. Dependence 
of critical strain on the specimen thickness shows that a criterion based on critical strain 
is not enough for predicting ductile fracture. Since the effect of thickness in fracture 
process can be accommodated through the crack tip triaxial Ity, it is better to have a 
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critit'ai «'untinuum parameter depending on both strain and triaxial ity to predict ductile 
crack growth initiation. The critical damage parameter [Dc), as identified in this thesis, 
ftiifills the above requirement. 




Thickness B = Smin 



° o Experimental 

Plane strain 

Three dimensional 



Load point deflection(niiii) 


Fig 6.2. Load deflection curve. 
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Fig 6.6. Variation of plastic strain along the 
thickness. 


2 4 6 8 10 12 14 16 

Thickness (mm) 

Fig 6.7. Fracture tougimess for specimens with 
different thickness. 





Fig 6.9. Critical strain for specimens with dif 
ferent thickness. 


Chapter 7 

SUMMARY AND CONCLUSIONS 


At microscopic level, a ductile fracture process depends on void nucleation. growth and 
coalescence. At continuum level, it depends on the critical values of several continuum 
parameters. The role of triaxiallity in void growth phenomenon had been identified long 
back. But it has not yet been established as a proper fracture parameter .either singly 
or in combination with other parameters. In this study an attempt has iDeen made to 
explain the void growth phenomenon in Mode I ductile fracture where triaxiallity plays 
a major role. AISI-1090, 1095 steels are chosen as the test materials which contain a 
large number of secondary particles. A damage growth law has been proposed with the 
help of Lemaitre’s continuum damage mechanics theory together with the experimental 
results of LeRoy et al. The damage variable D, quantifying the void growth, is physically 
interpreted as the area void fraction at a point in a plane. Thomason s limit load model 
for void coalescence is used for predicting micro crack initiation. 

The results of Chapter 4 show that the ductile failure depends on both the plastic 
strain and the triaxiality ( Fig 4.12 ). Since, the damage variable D is & combination of 
both the plastic strain and triaxiality, its critical value [Dc) is taken as the continuum 
parameter for predicting micro crack initiation. It is clear from Table 4.4 that the critical 
damage Dc is almost independent of the specimen geometry. Hence, Dc is taken as a 
material property. For AISI-1090 spherodised steel, it is coming out to be equal to 0.05. 
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(, !japt<‘r 5 is «T.n application of the critiral riamao-o 

damage as the continuum parameter 

and austenite grain size as the length parameter in a local criterion ( equation 5.1 ) for 
pred.ctmg the crtu-lt growth initiation in ductile materials. The numericallv predicted and 

the experintentally obtained values of fracture toughness ( J„) are in close agreement 

(s«* St*ction 5.5 ). 


The fractographs ( Figs .5.6(b), (c) ) of the fractured surface reveal that the critical 
stretch zone width I ) lies within the range of the measured grain size ( 30, .m - 

60, tm ). The computed value of { SZW ), is 32pm which is almost equal to the average 
austenite grain size ( IBpm ) of the material. The critical value of crack tip opening 
displacement ( [CTOD]^ } is obtained as 56pm from the computational results ( Fig 

0.17 ). Thus, {CTOD)g is 1.8 times {SZW)^ which seems reasonable. For a .semicircular 
crack tip, this value is likely to be 2 [61]. 

The J verses CTOD curve ( Fig 5.17 ) gives the slope of the curve as 'lAco- This 
result falls in the range of experimental results [ 57 ] for plane strain condition. 

The ratio of the plastic zone to damage zone ( Figs 5.13. 5.14 ) is obtained as 5.6. It 
shows damage as a highly localised event. Sometimes a fracture process zone is defined 
within the plastic zone and is of the size of the critical crack tip opening displacement 
[ 4 ]. A critical damage zone of one austenite grain size falls within the fracture process 
zone. 


As aaa overall conclusion of the results of Chapter 5, it can be stated that the 
proposed criterion ( equation 5.1 ) is reasonably good in predicting crack growth initiation 
in ductile materials. 

The results of Chapter 6 ( Figs 6.3 - 6.6 ) give the variations of continuum parameters 
at the crack tip along the thickness of the specimen. It also shows the dependence of 
critical value of Jj and critical plastic strain {elg)c on the thickness of the specimen ( Figs 
6.7 - 6.9 ). The qualitative trends of the results are consistent with the earlier published 
results [ 58 ]. 
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la the Ijgiil of liie above discussion the following conclusions can be drawn: 

1. The void growth law, proposed in this study for the class of materials like AISI-1090, 
AISM095 spherodised steel, is a reasonably good approximation. 

2. rhonuison s limit load model is a good representation of micro crack initiation due to 
void growth an<l audesccnce in ductile materials. 

3. Ductile fracture process is influenced by both the plastic strain and the triaxial Ity. 

4. The critical damage parameter (Dc) is a good measure of ductile failure initiation due 
to void growth. 

5. Choice of austenite grain size as the characteristic length parameter in crack growth 
initiation criterion leads to a good estimation of fracture parameters. 

6. The elastic-plastic FE analysis together with CDM can estimate fracture parameters 
reasonably well. 

Scope of Further Work 

The CDM theory together with FE analysis has been applied here to explain Mode 
I ductile fracture. The proposed damage growth law and the fracture criterion can be 
extended to other problems of fracture mechanics and metal forming processes, some 
important aspects are given below: 

1. Application of the proposed damage growth law and void coalescence condition in 
metal forming processes to predict the central burst event. 

2. Application of the proposed crack growth initiation criterion to mixed mode analysis 
of the ductile fracture. 
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APPENDIX 


Simulation of Necking 

Necking is a piu'noinenon which occurs due to localization of plastic deformation in a 
small region, hor example, in tension test of a cylindrical bar, a large reduction of 
diameter takes place in a small region due to local plastic flow where as the other part of 
the bar remains elastic. Definitely , it is an instability problem which can be treated in 
different ways. The common method to simulate instability is either to introduce a small 
perturbation in geometry or load ( imperfect approach [ 53 ] ) or to superpose on the 
di.splacement field of the critical load a part of the eigenmode ( perfect approach ). Here 
the imperfect approach of perturbing the geometry is used. It is done by introducing a 
small imperfection of the order of 0.1% in the diameter. 


The phenomenon is explained below for a one dimensional problem with the help of 
two elements ( Fig A1 ). The imperfection has been introduced at the bottom of element 
2, i.e. at node 3. For simplicity we assume that = 0 and Atd is the prescribed 
displacement at the end of the bar. One dimensional strain increments in elements (1) 
and (2) are 

An^ - Au^ (Al) 


Ae^ = 




and 


Ac' 


Au^ - Au^ Au^ 


since Au^ is zero. 


(A2) 


I2 h 

Here h and h are the lengths of the elements as shown in Fig Al(a). Since the element 
(2) has a lesser cross section than that of element (1), the stress level in element (2) is 
higher than that of element (1) to maintain equilibrium. These stress levels are indicated 
by points 1 and 2 on one dimensional stress-strain curve in Fig Al(b). Since the slope 
of the stress strain curve at point 2 is less than that at point 1, element (2) requires 
a greater strain increment than that of element (1) to have the same stress increment. 
More over, element (2) requires a larger stress increment to equilibrate the load since its 
cross section is less. Thus, for equilibrium it is required that 

Ac' > Ae^. (A3) 
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Since A«* is a fixe<i lirescrilnxi displacement increment, a situation will arise when 

(A4) 

I hen t'quiition I, AI) shows that element(l) experiences negative strain increment or un- 
loa<ling where ;ls the positive strain increment in element (2) increases. 
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